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ABSTRACT 



We present the results of a series of radio, optical. X-ray and y-ray observations of the BL Lac object S50716+714 carried out between 
April 2007 and January 2011. The multi-frequency observations were obtained using several ground and space based facilities. The 
' intense optical monitoring of the source reveals faster repetitive variations superimposed on a long-term variability trend at a time 

CO ' scale of ~ 350 days. Episodes of fast variability recur on time scales of ~ 60 - 70 days. The intense and simultaneous activity at 

optical and -y-ray frequencies favors the SSC mechanism for the production of the high-energy emission. Two major low-peaking 
radio flares were observed during this high optical/y-ray activity period. The radio flares are characterized by a rising and a decaying 
stage and are in agreement with the formation of a shock and its evolution. We found that the evolution of the radio flares requires 
' a geometrical variation in addition to intrinsic variations of the source. Different estimates yield a robust and self-consistent lower 

j_j I limits of (5 > 20 and equipartition magnetic field 6^, > 0.36 G. Causality arguments constrain the size of emission region 6 < 0.004 

, mas. We found a significant correlation between flux variations at radio frequencies with those at optical and y-rays. The optical/GeV 

■ ■ ■ flux variations lead the radio variability by ~ 65 days. The longer time delays between low-peaking radio outbursts and optical flares 

imply that optical flares are the precursors of radio ones. An orphan X-ray flare challenges the simple, one-zone emission models, 
rendering them too simple. Here we also describe the spectral energy distribution modeling of the source from simultaneous data 
taken through different activity periods. 

Key words, galaxies: active - BL Lacertae objects: individual: S5 0716+714 - Gamma rays: galaxies - X-rays: galaxies - radio 
continuum: galaxies 

1. Introduction 



* Member of the International Max Planck Research School (IMPRS) The BL Lac object S5 0716+714 is among the most extremely 
for Astronomy and Astrophysics at the Universities of Bonn and variable blazars. The optical continuum of the sot irce is so fea- 
Cologne tureless that it is hard to estimate its redshift. iNilsson et al.l 
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(I2008h claimed a value of z = 0.3 1 + 0.08 based on the photo- 
metric detection of the host galaxy. Very recently, the detection 
of intervening Lya systems in the ultra-violet spectrum of the 
source confir ms the earUer est i mates with a redshift value 0.23 15 
<z< 0.3407 (lDanforthetal.Ll20T2l) . This source has been clas- 
sified as an intermediate-peaked blazar (IBL) by iGiommi et al.l 
(Il999l) . as the frequency of the first spectral energy distribution 
(SED) peak varies between lO''* and 10'^ Hz, and thus does not 
fall into the wavebands specified by the usual definitions of low 
and high energy peaked blazars (i.e. LBLs and HBLs). 

S5 0716+714 is one of the brightest BL Lacs in the optical 
bands and has an optical intraday va riability (IDV) duty cycle of 
nearly one dWagner & Witzellll995h . Unsurprisingly, this source 
has been the subject of several optical monitoring campaigns on 



intraday (ID V) timescales (e.g. Wagner et al.L 1996t iRani et al.L 
[201 1: Monta gni et al.Ll2006HGupta et alil2009l" 1201 2L and refer- 



ences therein). The source has shown five major optical ou tbursts 
separated roughly by ~3.0+0.3 years (iRaiteri et al.Ll2003l) . High 
optical polarization ~ 20% - 29% has also b een observed in the 
source dTakalo et al l [l^lF^rretan.[T997h . lGuptaet al.l(l2009h 
analyzed the excellent intraday opti cal light curves of the source 
observed bv iMontagni et al.l (l2006h and reported good evidence 
of nearly periodic oscillations ranging between 25 and 73 min- 
utes on several diff'erent nights. Good evidence for the presence 
of ~ 15-mi n periodic oscillatio ns at optical frequencies has been 
reported bv lRani et al. IdMoi). A detailed multiband short-term 
o ptical flux and colo r variability study of the source is reported 
in lRani et al. l (l20T0bl) . There we found that the optical spectra of 
the source tend to be bluer with increasing brightness. 

There is a series of papers covering flux va riability studies 
dOuirrenbach et al.L 119911: IWagner et al.L I1996L and references 
there in) and morphologi c al/kinematic s tudies at radio frequen 



'1988^. 'Jorstad et al.L _ 

[19861 iBach et a'nT 2005 Rastorgueva et al.L 1201 iL and references 



cies dWitzel et a 



200 iL lAntonucci et all 



therein). Intraday variability at radio wavelengths is likely to be 
a mixture of intrinsic and extrinsic (due to interstellar scintil- 
lation) mechanisms. A significant correlation between optical - 
radio flux variations at day-to-day timescales has been reported 
by^agner et al. (1996). The frequency dependence of the vari- 
ability index at radio b ands is not found to be consistent with 
interstellar scintillation dFuhrmann et al.L l2008l) . which implies 
the presence of very small emitting regions within the source. 
However, the IDV time scale does show evidence in favor of 
an annual modulation, suggesting that the IDV of S5 0716+71 4 
could be dominated by interstellar scintillation dLiu et al 1 [MI- 
EGRET on board the Compton Gamma-ray Observatory 
( CGRO) detected high-energy y-ray (>100 MeV ) emission from 
S5 0716+714 several t imes from 1991 to 1996 dHartman et al.L 
[T999H L in et al I Il995h . Two strong y-ray flares on September 
and Octobe r 2007 were detected in the source with AGILE 
dChen et al.L l2008l) . These authors have also carried out 
SED modeling of the source with two synchrotron self- 
Compton (SSC) emitting components, representative of a slowly 
and a rapidly va riable componen t , respe ctively. Observations 
by BeppoSAX dTagliaferri et al.L l2003l) and XMM-Newton 
dFoschini et al.L l2006l) provide evidence for a concave X-ray 
spectrum in the 0.1 - 10 keV band, a signature of the pres- 
ence of both the steep tail of the synchrotron emission and the 
flat part of the Inverse Compton (IC) spectrum. Recently, the 
MAGIC collaboration reported the first detectio n of VHE y-rays 
from the source at the 5.8cr significance level (lA nderhub et al.L 
120091) . The discovery of S5 0716+714 as a VHE y-ray BL 
Lac object was triggered by its very high optical state, sug- 
gesting a possible correlation between VHE y-ray and the op- 



tical emissions. This source is also among the bright blazars 
in the Fermi/LAI (Large Area Telescope) Bright AGN Sample 
(LBAS) (Abdoet al., 2010a), whose GeV spectra are governed 
by a broken power law. The combined GeV - TeV spectrum of 
the source displays absorption-li ke features in 10 - 100 GeV en- 
ergy range dSenturk et al.L 1201 lb . 

The broadband flaring behavior of the source is even more 
complex. A literature study reveals that the broadband flar- 
ing activity of the source is not simultaneous at all frequen- 



cies dChen et a n.l2008HVillata et al.Ll2008HVittorini et a"n.l2009L 
lOstorero et al.L l2006l) . Also, it is hard to explain both the slow 
modes of variability at radio and hard X-ray bands and the 
rapid variability observed in the optical, soft X-r ay, and y-ray 



bands using a single component SSC model (seelVillata et al 



2008t IGiommi et all I2008L IChen et al.L I2008L IVittorini et al 



20091) . The X-ray spectrum of the source con tains contribu- 



tions from both synchrotr on and IC emission ([Foschini et al.L 
i2006l : iFerrero et al.L l2006l) and the simultan eous optical-GeV 
variations favor an S SC emission mechanism dChen et al.Ll2008L 
IVittorini et aLll2009l) . 

Despite of several eff'orts to understand the broadband flar- 
ing activity of the source, we still do not have a clear knowledge 
of the emission mechanisms responsible for its origin. In partic- 
ular, the location and the mechanism responsible for the high- 
energy emission and the relation between the variability at dif- 
ferent wavelengths are not yet well understood. Therefore, it is 
important to investigate whether a correlation exists between op- 
tical and radio emissions, which are both ascribed to synchrotron 
radiation from relativistic electrons in a plasma jet. If the same 
photons are up-scattered to high energies, simultaneous variabil- 
ity features could be expected at optical - GeV frequencies. But 
the observed variability often challenges such scenarios. To ex- 
plore the broadband variability features and to understand the 
underlying emission mechanism, an investigation of long-term 
variability over several decades of frequencies is crucial. The 
aim of the broadband variability study reported in this paper is 
to provide a general physical scenario which allows us to put 
each observed variation of the source across several decades of 
frequencies in a coherent context. 

Here, we report on a broadband variability study of the 
source spanning a time period of April 2007 to January 2011. 
The multi-frequency observations comprise GeV monitoring by 
Fermi/LALT, X-ray observations by Swift-XRT, as well as opti- 
cal and radio monitoring by several ground based telescopes. 
More explicitly, we investigate the correlation of y-ray activ- 
ity with the emission at lower frequencies, focusing on the indi- 
vidual flares observed between August 2008 and January 201 1. 
The evolution of radio (cm and mm) spectra is tested in the 
context of a standard shock-in-jet model. The broadband SED 
of the source is investigated using a one-zone synchrotron self- 
Compton (SSC) model and also with a hybrid SSC plus external 
radiation Compton model. In short, this study allows us to shed 
light on the broadband radio-to-y-ray flux and spectral variabil- 
ity during a flaring activity phase of the source over this period. 

This paper is structured as follows. Section 2 provides a brief 
description of the multi-frequency data we employed. In Section 
3, we report the statistical analysis and its results. Our discussion 
is given in Section 4 and we present our conclusions in Section 
5. 



2. Multi-frequency Data 

From April 2007 to February 20 11, S5 0716+714 was observed 
using both ground and space based observing facilities. These 
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Table 1. Ground based radio observatories 



Observatory 


Tel. dia. 


Frequency (GHz) 


Effelsberg, Germany 


100 m 


2.7, 4.8, 6.7 






8.3, 10.7, 15, 23 






32, 43 


UMKAU, Uj>A 


26 m 


4.0, 0, 14.3 


NOTO, Italy 


32 m 


5, 8, 22, 43 


Urumqi, China 


25 m 


4.8 


OVRO, USA 


40 m 


15 


Metsahovi, Finland 


14 m 


37 


PdBI, France 


6x15 m 


86, 143, 230 


Pico Veleta, Spain 


30 m 


86, 143, 230 


SMA, USA 


8x6 m 


230, 345 



multi-frequency observations of the source extend over a fre- 
quency range between radio and y-rays including optical and 
X-rays. In the following subsections, we summarize the obser- 
vations and data reduction. 



2.1. Radio and mm data 

We collected 2.7 to 230 GHz radio wavelength data of the source 
over a time period of April 2007 to January 201 1 (JD'Q = 200 
to 1600) using the 9 radio telescopes listed in Table 1. The 
cm/mm radio light curves of the source were observed as a 
part of observations within the framework of F-GAMMA pro- 
gram (Fermi-GST rela t ed monitoring p r ogram of y-ray blazars, 
iFuhrmann et~an 120071: lAngelakis etatl l2008l) . The overall fre- 
quency range spans from 2.7 GHz to 230 GHz using the 
Effelsberg 100 m telescope (2.7 to 43 GHz) and the IRAM 30 m 
Telescope at the Pico Veleta (PV) Observatory (86 to 230 GHz). 
These flux measurements were performed quasi-simultaneously 
using the cross-scan method slewing over the source position, 
in azimuth and elevation direction with adaptive number of sub- 
scans for reaching the desired sensitivity. Subsequently, atmo- 
spheric opacity correction, pointing off-set correction, gain cor- 
rection and sensitivity connection have been applied to the data. 

This source is also a part of an ongoing blazar mon- 
itoring program at 15 GHz at the Owens Valley Radio 
Observatory (OVRO) 40-m radio telescope which provides the 
radio data sampled at 15 GHz. We have also used the com- 
bined data from the University of Michigan Radio Astronomy 
Observatory (UMRAO; 4.8, 8 and 14.5 GHz. lAller et"anil985l) 
and the Metsahovi Radio Ob servatory (MRO; 22 and 37 GHz; 
iTerasranta et al.L 1 19981 l2004l) . which provide us with radio light 
curves at 5, 8, 15 and 37 GHz. Additional flux monitoring at 
5, 8, 22 and 43 GHz radio bands is obtained using 32 m tele- 
scope at NOTO radio observatory. The Urumqi 25 m radio tele- 
scope monitors the source at 5 GHz. The 230 and 345 GHz 
data are provided by the Submillimete r Array (SMA) Observer 
Centefl data base (iGurwell et al.Ll2007h . complemented by some 
measurements from PV and the Plateau de Bure Interferometer 
(PdBI). The radio light curves of the source are shown in Fig 
[U The mm observations are closely coordinated with the more 
general flux density monitoring conducted by IRAM, and data 
from both programs are also included. 



' JD' = JD-2454000 

^ http://smal .sma.hawaii.edu/callist/callist.html 



2.2. Optical data 

Optical V passband data of the source were obtained from the 
observations at the 1.5-m Kanata Telescope located on Higashi- 
Hiroshima Observatory over a time period of February 14, 
2009 to June 01, 2010 (JD ' = 877 to 1349) The Triple Range 
Imager and SPECtrograph (IWatanabe et al.L l2005l) was used for 
the observations from JD' = 612 to 1228. Th en, the HOWPol 
(Hiroshima One-shot Wide-field Polarimeter, iKawabata et all 
|2008) was used from JD' 1434 to 2233 with the V-band filter. 
Exposure times for an image ranged from 10 to 80 s, depending 
on the magnitude of the object and the condition of sky. The pho- 
tometry on the CCD images was performed in a standard proce- 
dure: after bias subtraction and flat-field division, the magnitudes 
were calculated using the aperture photometry technique. 

Additional optical V-passband data were obtained from the 
2.3 m Bok Telescope of Steward Observatory from April 28, 
2009 through June 2, 2010 (JD' = 950-1350). These data are 
from the public data archive that provides the results of po- 
larization and flux monitoring of bright y-ray blazars selected 
from the ferm//LAT-monitored blazar lisl0. We have also in- 
cluded optical V passband archival data extracted from the 
American Association of Variable Star Observers (AAVSO; 
see http://www.aavso.org/ for more information) over a period 
September 2008 to January 2011 (JD' = 710-1600). The com- 
bined optical V passband flux light curve of S5 0716-1-714 is 
shown in Figure |2](c). 

2.3. X-ray data 

The X-ray (0.3 - 10 keV) data were obtained by Swift-XSJ 
over a time period of September 2008 to January 2011 (JD' = 
710-1600). The Swift-XKl data were processed using the most 
recent versions of the standard Swift tools: Swift Software ver- 
sio n 3.8, FTOOLS v ersion 6.11 and XSPEC version 12.7.0 (see 
iKalberla et al.Ll2005L and references therein). 

All of the observations were obtained in photon counting 
(PC) mode. Circular and annular regions are used to describe 
the source and background areas, respectively, and the radii 
of both regions depend on the measured count rate using the 
FTOOLS script xrtgrblc. Spectral fitting was done with an ab- 
sorbed power-l aw with Nh - . 31 x l O^'cm^^ set to the Galactic 
value found bv IKalberla et al.1 (l2005h . For three of the observa- 
tions, there were more than 100 counts in the source region, and 
a chi-squared statistic is used with a minimum bin size of 20 
cts/bin. For the bin centered on JD' = 1214, only 62 counts were 
found in the source region, and the unbinne d data are fitted us- 
ing C-statistics as described bv lCashl([T979h . One sigma errors 
in the de-absorbed flux were calculated assuming that they share 
the same percentage errors as the absorbed flux for the same time 
and energy range. The X-ray light curve of the source is shown 
inFig.Hb). 

2.4. y-ray data 

We employ y-ray (100 MeV - 300 GeV) data collected be- 
tween a time period August 08, 2008 to January 30, 2011 
(JD'= 686-1592) in survey mode by FermZ/LAT. The LAT data 
are analyzed using the standard ScienceTools (software version 
v9.23.1) and the instrument response function P7V(fl Photons 
in the Source event class are selected for this analysis. We se- 



^ http://james.as.arizona.edu/~psmith/Fermi 
http://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/overview.html 
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Fig. 1. Radio to mm wavelength light curves of S5 0716+714 observed over the past ~3 years. For clarity, the light curves at different 
frequencies are shown with arbitrary offsets (indicated by a "Frequency + x Jy" label). The major radio flares are labeled as "RO", 
"R6" and "R8" (see Fig.|2]for the details of labeling). 



lect y-rays with zenith angles less than 100 deg to avoid con- 
tamination from y-rays produced by cosmic ray interactions in 
the upper atmosphere. The diffuse emission from our Galaxy is 
modeled using a spatial modeQ (gal_2yearp7v6_v0.fits) which 
was derived with Fermi/LAT data taken during the first two 
years of operation. The extragalactic diffuse and residual in- 
strumental backgrounds are modeled as an isotropic component 



(isotropic_p7v6source.txt) with both flux and spectral photon in- 
dex left free in the model fit. The data ana lysis is done with an 
unbinned maximum likelihood technique (Ma ttox et al.L 1 19961) 
using the publicly available toolfl We analyzed a Region of 
Interest (Rol) of 10° in radius, centered on the position of the 
y-ray source associated with S5 0716+714, using the maximum- 
likelihood algorithm implemented in gtUke. In the Rol model. 



http://fermi.gsfc.nasa.gov/ssc/data/access/lat/BackgroundModels.html * http://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/likelihood_tutorial.html 
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Fig. 2. Light curves of 0716+714 from y-ray to radio wavelengths (a): GeV light curve at E>100 MeV, (b): X-ray light curve at 
0.3-10 KeV, (c): optical V passband light curve and (d) : 5 to 230 GHz radio light curves. Vertical dotted lines are marked w.r.t. 
different optical flares labeling the broadband flares as "G" for y-rays, "X" for X-rays, "O" for optical and "R" for Radio followed 
by the number close to flare. The yellow area represents the period for which we construct the broadband SEDs of the source (see 
Section [33] for details). 



we include all the 24 sources within 10° with their model pa- 
rameters fixed to their 2FGL catalog values, as none o f these 
sources are reported to be variable (see lNolan et al.[|2012[ for de- 
tails). The contribution of these 24 sources within the Rol model 
to the observed variability of the source is negligible as they 
are very faint compared to S5 0716+714. The LAT instrument- 
induced variability is tested with bright pulsars and is ~2% and 



is much smaller than the statistical errors reported for the source 
dAckermann et al.Ll20T2h . 

Source variability is investigated by producing the light 
curves (E > 100 MeV) by likelihood analysis with time bins of 
1, 7 and 30 days. The bin-to-bin exposure time variation is less 
than 7%. The light curves are produced by modeling the spec- 
tra over each bin by a simple power law which provide a good fit 
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Table 2. Variability time scales at radio wavelengths 



Frequency 




t,,,,-! (days) 


pi 


t,o,-2 (days) 


15 GHz 


0.95+0.03 


100+5 


1.32+0.04 


195+5 


37 GHz 


1.50+0.13 


100+5 


1.64+0.10 


200+5 


43 GHz 


0.99+0.02 


90+5 


1.23+0.04 


180+5 


86 GHz 


1.60+0.10 


90+5 


0.89+0.02 


180+5 


230 GHz 


1.04+0.03 


90+5 


1.31+0.03 


180+5 



* : P(f) ~ /'', p is the slope of the power law fit. 



over these small bins of time. Here, we set both the photon index 
and integral flux as free parameters. We will use the weekly av- 
eraged light curves for the multi-frequency analysis, as the daily 
averaged flux points have a low TS (Test Statistics) value (<9). 
In a similar way, we construct the GeV spectrum of the source 
over different time periods (see Section l33] for details). We split 
the 0.1 to 100 GeV spectra into 5 different energy bins : 0.1 - 
0.3, 0.3 - 1.0, 1.0 - 3.0, 3.0 - 10.0 and 10.0 - 100.0 GeV. A sim- 
ple power law with constant photon index F (the best fit value 
obtained for the entire energy range) provides a good estimate 
of integral flux over each energy bin. The GeV spectrum of the 
source is investigated over four different time periods, which rep- 
resent different brightness states of the source and will be used 
in broadband spectral modeling. The details of the broadband 
spectral analysis are given in Section l33l 

3. Analysis and Results 

3. 1 . Light Curve Analysis 
3.1.1. Radio frequencies 

Fig. [Udisplays the 2.7 - 230 GHz radio band light curves of the 
source observed over the past ~ 3 years. The source exhibits sig- 
nificant variability, being more rapid and pronounced towards 
higher frequencies over this period with two major outbursts. 
Towards lower frequencies (<10 GHz) the individual flares ap- 
pear less pronounced and broaden in time. 

To quantify the strength of variability at different radio fre- 
quencies, we extract the time scale of variability {t^ir) from the 
observed light curves. W e employed the structure fu nction (SF) 
( Simonetti et al.L Il985h analysis method following iRani et al] 
(I2009h . The radio SF curves are shown in Fig. |3] 

At 15 GHz and higher radio frequencies, the SF curve fol- 
low a continuous rising trend showing a peak at t^ari, following 
a plateau and again reaching a maximum at f,,„r2. However, the 
SF curve at 10 GHz and lower frequencies do not reveal a sharp 
break in the slope as the variability features seem to be milled 
out at these frequencies. We do not consider variability features 
at time lags longer than half of the length of the observations due 
to the increasing statistical uncertainty of the SF values in this re- 
gion. To extract the variability timescales, we fitted a power law 
(P(/) ~ /^) to the two rising parts of the SF curves. The variabil- 
ity timescale is given by a break in the slope of the power-law 
fits. In Fig.[3](top), the red curves represent the fitted power-law 
to rising trend of SF curves and the vertical dotted lines stand for 
tvari and f,,„r2. The best fitted values of the power-law slopes and 
the estimated timescale of variability are given in Table |2] Thus, 
the SF curve reveals two different variability time scales, one 
which reflects the short-term variability (fvari) while the other 
one refers to the long-term variability (tyari)- 

The SF value is proportional to the square of the amplitude 
of variability, so to compare the strength of the variability at 
different frequencies, we produce SF vs frequency plots at 
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Time lag (days) 




10 100 1000 



Time lag (days) 

Fig. 3. Top: Structure functions at radio frequencies. The solid 
curves represent the best fitted power laws. The dotted lines 
in each plot indicate the timescale of variability (?,«,). Bottom: 
Structure functions at optical and y-ray frequencies. 

different time lags. In fig. |4l we show the SF vs frequency plot 
at a time lag of 100 days. The source displays more pronounced 
and faster variability at higher radio frequencies compared 
to lower ones, peaking at a frequency of 43 GHz, and have 
similar amplitude at higher frequencies. It seems that the radio 
variability saturates above this frequency. We notice a similar 
trend at 50 and 200 day time lags. Thus, for different time lags 
the variability strength exhibits a similar frequency dependence. 



3.1.2. Optical frequencies 

The source exhibits multi-flaring behavior at optical frequencies 
with each flare roughly separated by 60 - 70 days (see Fig.|2jc)). 
The optical V band SF curve in Fig. [3] (bottom) shows rapid 
variability with multiple cycles of rises and declines. The first 
peak in SF curve appears at a timescale ~ 30 days which 
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Fig. 4. SF vs frequency at GHz frequencies for a structure func- 
tion time = 100 days. 
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Fig. 5. LSP analysis curve showing a peak at a period of 63 and 
359 days. 



is followed by a dip at ~ 60 days. This peak corresponds to the 
fastest variability timescale. The peaks in the optical SF curve at 
tyar - 90, 180, 240 and 300 days represent the long-term vari- 
ability timescales. This indicates a possible superposition of a 
short 30-40 day time scale variability with the long-term vari- 
ability trend. 

Multiple cycles in the optical SF curve represent the nearly 
periodic variations at ~60 da ys time s cale. We used the Lomb- 
Scargle Periodogram (LSP) dLombl [T976t iRani et al.L l2009l) 
method to test the presence of this harmonic. The LSP analy- 
sis of the whole data set is displayed in Fig.|5] The LSP analysis 
reveals two significant (>99.9999% confidence level) peaks at 
359 and 63 days. The peak at 359 days is close to half of the 
duration of observations, so it is hard to claim this frequency 
due to limited number of cycles. The nearly annual periodicity 
could also be effected by systematic effects due to annual ob- 
serving cycle. A visual inspection of the light curve indicates a 
total of 7 rapid flares separated by 60 - 70 days. Also, the LSP is 
only sensitive for a dominant white-noise process {Pj^if) oc /°). 
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Fig. 6. Optical V passband Hght curve of S5 0716+714 with dif- 
ferent time binning. The light curve appears as a superposition 
of fast flares on a modulated base level varying on a (350+9) day 
timescale. These slower variations can be clearly seen in 75 days 
binned light curve (error bar represents variations in flux o ver the 
binned period). The dashed line represents a spline interpolation 
through the 75 day binned light curve. Dotted lines are obtained 
by shifting the spline by +0.65 mag. 



It is for this reason that we further inspect the significance of 
this freq uency with the P ower Spectral Density (PSD) analysis 
method ( IVaughanl l2005h which is a powerful tool to search for 
periodic signals in time series, including those contaminated by 
white noise and/or red noise. To achieve a uniform sampling in 
the optical data, we adopt a time-average binning of 3 days. We 
found that the significance of the period at ~ 60 days is only 2 
cr. We therefore can not claim a significant detection of a quasi- 
periodic variability feature at this frequency. 

We also notice that during the course of our optical monitor- 
ing, the peak-to-peak amplitude of the short-term variations re- 
mains almost constant, ~ 1.3 magnitudes. Hence, the variability 
trend traced by the minima is very similar to that by the maxima 
of the light curve during the course of ~ 2 years of our moni- 
toring. The constant variability trend is displayed in Fig. |6] In 
this figure, the dashed line denotes a spline through the 75-days 
binned light curve and the dotted lines are obtained by shifting 
the spline by +0.65 mag. So, the observed variations fall within a 
constant variation area. A constant variability amplitude in mag- 
nitudes implies that the flux variation amplitude is proportional 
to the flux level itself (following m\ - ml oc logn){f\ j f2)). This 
can be easily interpreted in terms of va riations of th e Dopp ler 
boosting factor, 6 = [r(l - p cos 6*)] ' dRaiteri et al.Ll2003l) . In 
such a scenario, the observed flux is relativistically boosted by 
a factor of 5^ and requires a variation in (5 by a factor of ~ 1.2. 
Such a change in 6 can be due to either a viewing angle {6) vari- 
ation or a change of the bulk Lorentz factor (F) or may be a 
combination of both. A change in 5 by a factor of 1 .2 can be 
easily interpreted as a few degree variation in 6 while it requires 
a noticeable change of the bulk Lorentz factor. Therefore, it is 
more likely that the long-term flux base-level modulations are 
dominated by a geometrical effect than by an energetic one. 

Hence, we consider that the optical variability amplitude 
remains almost constant during our observation s. A similar 
variability trend was also observed in this source bv lRaiteri et al.l 
( l2003h . They also noticed a possibility of nearly periodic oscil- 
lations at a timescale of 3.0 + 0.3 years, but due to the limited 
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time coverage this remains uncertain. The optical light curve of 
the source also displays fast flares with a rising timescales of 
~ 30 day. However, the rising timescale of the radio flares is of 
the order of ~ 60 days (see Table|5]l. Thus, the optical variability 
features are very rapid compared to those at radio wavelengths. 



3.1 .3. X-ray frequencies 

Fig.lab) displays the 0.3 — 10 keV light curve of S5 0716+714. 
Although the X-ray light curve is not as well sampled as those 
at other frequencies, the data indicate a flare at keV energies 
between JD' - 1 122 to 1 165. Due to large gap in the data train, 
it is not possible to locate the exact peak time of this flare. 
Still, it is interesting to note that this orphan X-ray flare is not 
simultaneous with any of the GeV flares nor with the optical 
flares. Its occurrence coincides with the optical - GeV minimum 
after the major flares at these frequencies (05/G5, see Fig.|2|. 



3.1.4. Gamma-ray frequencies 

The GeV light curve observed by FermifLP^ is extracted over 
the time period of August 2008 to February 201 1. Fig. Q shows 
the weekly and monthly averaged y-ray light curves integrated 
over the energy range 100 MeV to 300 GeV. There is a signif- 
icant enhancement in the weekly averaged y-ray flux over the 
time period of JD' = 900 to 1 1 10, peaking at JD' ~ 1 1 10, with 
a peak flux value of (0.57 + 0.05) x 10"* ph cm^^s"', which is 
~ 6 times brighter than the minimum flux and ~ 3 times brighter 
than the average flux value. Later it decays, reaching a minimum 
at JD' = 1 150, and then it remains in a quiescent state until JD' 
= 1200. The quiescent state is followed by another sequence of 
flares. The source displays similar variability features in the con- 
stant uncertainty light curve obtained using th e adaptive binning 
analysis method following iLott et all (|2012|) . A more detailed 
discussion of GeV variability is given in Rani et al. (2012). 

The source exhibits a marginal spectral softening during the 
quiescent state in the monthly averaged light curve. The y-ray 
photon index (F) changes from ~ 2.08 + 0.03 (II) to 2.16 + 0.02 
(III), then again to 2.04 + 0.04 (IV). Diff'erent activity states of 
the source are separated by vertical lines in Fig. |7] We notice no 
clear spectral variation in the weekly averaged light curve due 
to large uncertainty and scatter of individual data points. 

The SF curve of the source at GeV frequencies is shown in 
Fig. [3] (bottom). The variability timescales are extracted using 
the power-law fitting method as described above, which gives a 
break at ~ 30 and 180 days. We notice that the variability fea- 
tures at tyar ~ 180 days are also observed at radio and optical 
frequencies. However, the faster variability (t^ar ~ 30 days) ob- 
served at optical and y-ray frequencies does not extend to radio 
wavelengths. A similar long-term variability timescale at y-ray, 
optical and radio frequencies provides a possible hint of a co- 
spatial origin. In the following sections, we will search for such 
possible correlations. 



3.2. Correlation Analysis 

The multi-frequency light curves of S5 0716+714 are presented 
in Figure |2] The analysis presented in this section is focused 
on a time period between JD' = 800 to 1400, which covers the 
two major radio flares and the respective optical-to-y-ray flaring 




700 800 900 1000 1100 1200 1300 1400 1500 1600 
JD [2454000 +] 

Fig. 7. Gamma-ray flux and photon index light curve of 
S5 0716+714 measured with the Fermi/LAT for a time period 
between August 04, 2008 to February 07, 2011. The blue sym- 
bols show the weekly averaged flux while monthly averaged val- 
ues are plotted in red. Diff'erent activity states of the source are 
separated by vertical green lines. A marginal softening of spec- 
trum over the quiescent state can be seen here. 



activity. The source displayed multiple flares across the whole 
electromagnetic spectrum over this period, which we label as 
follows. We mark the vertical dotted lines w.rt. different optical 
flares, labeling the broadband flares as "G" for y-rays, "X" for X- 
rays, "O" for optical and "R" for radio followed by the number 
adjacent to them. For example, the optical flares should be read 
as "Ol" to "09". 

To search for possible time lags and to quantify the correla- 
tion among the multi-frequency light curves of the source, we 
computed discrete cross-c orrelation functions (DCF s) following 
the method described by (lEdelson & Krolikl Il988h. The details 
of this method are also discussed bv iRani et al.l ( |2009|) . In the 
following sections, we will discuss in detail the correlation be- 
tween the observed light curves. 



3.2.1. Radio correlation 

At radio wavelengths, the source exhibits significant flux vari- 
ability, being more rapid and pronounced towards higher fre- 
quencies. We label the two major outbursts as "R6" and "R8" 
(see Fig. |2]i. Apparently, the mm flares are observed a few days 
earlier than the cm flares. Our dense frequency coverage facili- 
tates a cross-correlation analysis between the different observing 
bands. Owing to its better time sampling, we have chosen the 
15 GHz light curve as a reference. Fig.|8]shows the DCF curves 
adopting a binning of 10 days. 

We used a Gaussian profile fitting technique to estimate the 
time lag and respective cross-correlation coefficient value for the 
DCF curves. Around the peak, the DCF curve as a function of 
time lag t can be reasonably well described by a Gaussian of 

the form: DCF{t) - ax exp[^^^j^], where a is the peak value 
of the DCF, b is the time lag at which the DCF peaks and c 
characterizes the width of the Gaussian function. The calculated 
parameter (a, b and c) values for each frequency are listed in 
Table |3] The solid curve represents the fitted Gaussian function 
in Fig. m 
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Fig. 8. The DCF curves among the different radio frequency light curves. The solid curves represent the best-fit Gaussian function. 



Table 3. Correlation analysis results among radio frequencies 



Frq. (GHz) 


a 


b (days) 


c (days) 


230 vs 15 


1.56i 


0.15 


7.96±2.23 


19.81+2.24 


86 vs 15 


0.944 


:0.13 


6.65±3.28 


25.80+4.28 


43 vs 15 


1.04i 


0.11 


5.95±2.08 


23.86+3.09 


37 vs 15 


1.13i 


0.09 


4.95+2.21 


29.39+2.81 


23 vs 15 


1.17i 


0.10 


3.74+1.50 


25.00+2.50 


10 vs 15 


0.89i 


0.09 


-1.01+1.09 


35.07+4.49 


8 vs 15 


0.84i 


0.08 


-1.09+1.01 


35.96+4.10 


5 vs 15 


0.84i 


0.10 


-1.23+1.25 


33.13+4.15 


2.72 vs 15 


0.59i 


0.12 


-78.75+12.39 


53.54+13.86 



a : peak value of the DCF, 

b : time lag at which the DCF peaks, and 

c : width of the Gaussian function (see text for details) 



Our analysis confirms the existence of a significant corre- 
lation across all observed radio-band light curves till 15 GHz 
with formal delays listed in Table |3] (parameter b). However, no 
pronounced delayed flux variations are observed at 10 GHz and 
lower frequencies. Also, the flux variations at 2.7 GHz seem to 



be less correlated with those at 15 GHz, which is obvious as the 
flaring behavior is not clearly visible below 15 GHz. 

The long term radio light curve shows three major radio 
flares, labeled as RO, R6 and R8 in Fig. [T] In the correlation 
analysis of the entire light curves, these flares are blended and 
folded into a single DCF. In order to separate the flares from 
each other, we performed a correlation analysis over three dif- 
ferent time bins which cover the time ranges of the individual 
radio flares: JD' = 500 to 750 (RO flare), JD' = 1000 to 1210 
(R6 flare) and JD' = 1210 to 1400 (R8 flare). The time lags of 
these flares relative to the 15 GHz data are estimated as above. 
However due to sparse data sampling, it was not possible to es- 
timate the time lags for R8 directly using the DCF method. So, 
for this flare, we first interpolate the data through a spline func- 
tion and then perform the DCF analysis, except at 23 GHz and 
33 GHz (due to long data gaps). The calculated time lags of each 
flare are given in Table |4] 

In Fig. |9l we report the calculated time lags as a function 
of frequency with 15 GHz as the reference frequency. The 
estimated time lag using the entire light curves are shown 
with blue circles. As we see in Fig. |9] the time lag increases 
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Table 4. Radio correlation analysis results for individual flares 



Frequency 




Time lag* (days) 




(GHz) 


RO Flare 


R6 Flare 


R8 Flare 


23 


2.9±1.4 


4.06±1.6 




33 


2.1±1.0 


5.88±2.1 




37 


4.2±2.6 


3.15±2.1 


4.0±1.0 


86 


4.2±1.0 


6.15±2.6 


5.0±1.8 


143 


5.8±1.8 


6.82±2.5 


7.0±1.5 


230 


6.0±2.4 


8.54±1.9 


9.0±2.0 



relative to the respective flares at 15 GHz 
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Fig. 10. The modeled radio flare, R6. The blue points are the 
observed data while the red curve represent the fit. 

Table 5. Fitted model parameters for R6 flare 



Frq. 




f,,,,,.. 


t,- 


t<; 


to 


GHz 








JD' [JD-2454000] 


15 


0.02±0.07 


4.15±0.14 


61.4±6.2 


37.9±3.5 


1191.4 ±0.9 


23 


0.71±0.23 


5.30±0.15 


32.3±4.8 


17.3±2.1 


1190.2 ±0.1 


37 


0.58±0.11 


8.20±0.59 


55.5±11.0 


18.6±3.2 


1189.1 ±0.7 


43 


0.45±0.35 


9.50±0.62 


60.5±9.4 


20.1±2.9 


1188.0 ±0.8 


86 


0.64±0.18 


10.6±2.48 


60.9±28.8 


25.1±25.6 


1186.0 ±0.5 


230 


0.72±0.11 


12.64±0.29 


50.3±2.4 


9.9 ±0.6 


1184.2 ±0.4 


* : see text for the extension of labels. 



Frequency (GHz) 



Fig. 9. The plot of time lag vs frequency, using 15 GHz as the 
reference frequency. Time lag vs frequency curves for individual 
flares are shown with different colors. The solid lines represent 
the best fitted power law in each case. 

with an increase in frequency and seems to follow a power 
law. Consequently, we fitted a power law, P{f) - Nf to 
the time lag vs frequency curve. The best fitted parameters 
are A? = 1.71 + 0.43, a = 0.30 + 0.08. For the individual 
flares, we find; = 1.07 + 0.06, a = 0.32 + 0.01 for the RO, 
N =\A5 + 0.61, a = 0.32 + 0.08 for R6, and N ^133+ 0.01, 
a = 0.29 + 0.03 for R8. We conclude that a common trend in the 
time lag (with 15 GHz as the reference frequency) vs. frequency 
curve is seen for all the three radio flares (RO, R6 and R8) with 
an average slope of 0.30. 

We also followed an alternative approach to estimate the time 
shift of the radio flares at each frequency w.rt. 15 GHz. The 
flares at each frequency can be modeled with an exponential rise 
and decay of the form: 

fit) = /o + fma, exp[(t - fo)/f,], for t < to, and (1) 

= /o + fnax exp[-(t - to)/ td] for t > to 

where fo is the background flux level that stays constant over the 
corresponding interval, f^ax is the amplitude of the flare, to is the 
epoch of the peak, and t, and td are the rise and decay time scales, 
respectively. Since R6 is the most pronounced and best sampled 
flare, we model this flare in order to cross-check the frequency 



vs time lags results obtained by the DCF method. As the flaring 
behavior is not clearly visible below 15 GHz, so we restrict this 
analysis to frequencies above 15 GHz. As there is no observation 
available during the flaring epoch at 23 and 86 GHz, we fix t^ and 
trf using the fit parameters from the adjacent frequencies. The 
best fit of the function f(t) for the R6 flare is shown in Fig. [10] 
and the parameters are given in Table |5] The estimated time shift 
around the R6 flare at each frequency w.rt. 15 GHz are shown in 
Fig.|9](red symbols) and the fitted power law parameters aie N = 
1.17 + 0.13, a - -0.31 + 0.03. Thus, this alternative estimate of 
the power law slope using the model fitting technique confirms 
the results obtained by the DCF analysis. 

3.2.2. Radio vs Optical correlation 

S5 0716+714 exhibits multiple flares at optical frequencies. The 
flares are roughly separated by 60 - 70 days. We labeled the dif- 
ferent optical flares as Ol - 09 as shown in Fig. |2] During this 
multi-flaring activity period two major flares are observed at ra- 
dio wavelengths. The radio flare R6 apparently coincides in time 
with 06 and R8 with 08. To investigate the possible correlation 
among the flux variations at optical and radio frequencies, we 
perform a DCF analysis using the 2-year-long simultaneous op- 
tical and radio data trains between JD' = 680 to 1600 (see Fig. 
|2]i. Note that the strength of flux variability increases towards 
higher frequencies, peaking at 43 GHz (see Fig. |3]l. Therefore, 
we choose two radio frequencies, 37 GHz and 230 GHz, in order 
to compare the strength of radio - optical correlation above and 
below the saturation frequency (43 GHz). The optical vs radio 
DCF analysis curves are shown in Fig.[TT](a). Multiple peaks in 
the DCF may reflect a QPO behavior at optical frequencies. As 
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Table 6. Fitted model parameters for R6 flare 



V vs 230 GHz 


V v,s 37 GHz 


lag (days) 


DCF Peak value 


lag (days) 


DCF Peak value 


-4±2.5 


0.43±0.10 


-2±2.5 


0.28+0.11 


63±2.5 


0.83±0.11 


66±2.5 


0.76±0.08 


120±2.5 


0.60±0.08 


124±2.5 


0.60±0.09 


181±2.5 


0.51 ±0.07 


183±2.5 


0.49±0.08 



the formal errors, we use half of the binning time. We summa- 
rize the optical vs 230 GHz and 37 GHz DCF analysis results in 
Table |6l 

The maximum correlation of the optical V passband with 
the 230 GHz light curve occurs at a 65 day time lag. However 
a second peak with lower peak coefficient also occurs close to 
zero time lag (see Fig. fTTT l. The analysis shows that the cross- 
correlation coefficient of the simultaneous radio - optical flare 
peaks 06-R6 and 08-R8 is lower than the cross-correlation co- 
efficient of the 05-R6 and 07-R8 flare peaks. In both cases the 
optical flares 05 and 07 are observed ~ 65 days earlier than 
the radio flares R6 and R8, respectively. The observed time lag 
between the optical and radio flares is consistent with the ex- 
trapolated frequency dependent time shift (as shown in Fig.|9]l to 
optical wavelengths. In Section l4n we will discuss this in detail. 

In order to quantify the correlation among optical and radio 
data, we generate flux - flux plots, which are shown in Fig. [TT] 
(b) - (c). For the following analysis we used a 1 day binning. 
Fig. [111(c) shows the time shifted 230 GHz (t-63) and 37 GHz 
(t-66) flux plotted vs the optical V-band flux. The time-shifted 
radio and optical V-band fluxes fall on a straight line, indicating 
a correlation. A Pearson correlation analysis reveals a significant 
correlation between the two data trains. We obtain the following 
values: rp = 0.59 and 99.93 % confidence level for 230 GHz 
(t-65) vs V-band and rp - 0.43 and 99.3 % confidence level for 
37 GHz (t-65) vs V-band, where rp is the linear Pearson correla- 
tion coefficient. Thus, we found a significant correlation among 
the time shifted radio vs optical V-band flux at a confidence level 
> 99 %. 

In contrast to this, the radio (with no time shift) vs optical 
V-band correlations are found to be not significant. Fig. [TT](c) 
shows 230 GHz and 37 GHz vs V-band flux-flux plots, and the 
correlation statistics are: 230 GHz vs V-band: rp = 0.40, 91 % 
confidence level and 37 GHz vs V-band: rp - 0.15, 74 % con- 
fidence level. Thus, the confidence level of the correlations is 
lower than 95 % in these cases. We also check the significance of 
the correlation statistics with a time shift of ~ 120 and 180 days 
and do not find a correlation to have a significance greater than 
2cr in any case. 

Hence, using DCF and linear Pearson correlation statistics, 
we have found a significant correlation among the flux variations 
at optical and radio frequencies with the optical V-band leading 
the radio fluxes at 230 and 37 GHz by ~ 63 and ~ 66 days, re- 
spectively. We therefore conclude that flux variations at optical 
and radio frequencies are correlated such that the optical vari- 
ability is leading the radio with a time lag of about two months. 

3.2.3. Radio vs Gamma-ray correlation 

We apply the DCF analysis method to investigate a possible cor- 
relation among flux variations at radio and y-ray frequencies. In 
Fig. [121 we report the DCF analysis results of the weekly av- 
eraged y-ray light curve with the 230 GHz radio data with a 
time bin size of 1 1 days. As the flux variations at 37 GHz are 
delayed by ~ 3 days w.rt. those at 230 GHz, we only show the 
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Fig. 11. (a) The DCF curve for optical V passband vs 37 GHz 
(in blue) and 230 GHz (in red) flux with a bin size of 5 days, 
(b) Radio flux vs optical V-band flux, (c) Time shifted radio flux 
vs optical V-band flux. The blue symbols show the time shifted 
230 GHz (t-65 days) data while 37 GHz (t-68 days) data are 
shown in red. 



DCF analysis curve w.rt 230 GHz. To estimate the possible peak 
DCF value and respective time lag, we fit a Gaussian function to 
the DCF curve with a bin size of 1 1 days. The best-fit function 
is shown in Fig. [T2land the fit parameters are a = 0.94 + 0.30, 
b - {67 ±3) days and c = (7 + 2) days. This indicates a clear cor- 
relation between the y-ray and 230 GHz radio light curves of the 
source with the GeV flare leading the radio flare by (67+3) days. 

To check the significance of the y-ray vs radio correlation, 
we produce flux-flux plots of the time shifted radio vs y-ray flux. 
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^ ' ' \ ' ' ] Table 7. Optical vs. y-ray cross-correlation analysis results 



40 80 120 

Time lag (days) 



3 




0.1 
Gamma 



0.3 0.3 0.4 0.5 

-ray Flux (le~' ph cm"^ 



Fig. 12. Top: DCF curve of the y-ray light curve w.r.t. the 
230 GHz radio light curve. The solid curve is the best fitted 
Gaussian function to the 11 day binned DCF curve. Bottom: 
Flux-flux plot of the shifted radio vs y-ray data. The blue sym- 
bols show the time shifted 230 GHz data while 37 GHz data are 
shown in red. 



Since the y-ray flux is weekly averaged, we use a time binning 
equal to seven days. The weekly averaged flux-flux plots of the 
time shifted 230 GHz (t + 67 days) and 37 GHz (t + 70 days) vs 
y-ray are shown in Fig. [12] (bottom) and the correlation statistics 
are: 230 GHz (t + 67 days) vs y-ray: rp = 0.37, 97.7 % confi- 
dence level and 37 GHz (t + 70 days) vs y-ray: rp = 0.33, 97.3 % 
confidence level. Thus, in each case the confidence level of the 
correlation is higher than 95 %. This supports a possible corre- 
lation among the flux variations at y-ray and radio frequencies 
with y-rays leading the radio emission by ~ 67 days. We also 
note that the time shifts are very similar to the time shifts ob- 
served between radio and optical bands (see Section l3.2.2) . We 
therefore expect a very short or no time delay between the flux 
variations at optical and y-ray frequencies. This will be investi- 
gated in the next section. 



Case 


Time duration 


Peak DCF value 


Time lag 




JD' [JD-2454000] 




days 


A 


total (840 - 1350) 


0.50±0.04 


0±5 


B 


removing 06 (1 150 - 1220) flare 


0.61 ±0.04 


1±5 


C 


before 06 flare (840 - 1150) 


0.80±0.08 


3±5 



3.2.4. Optical vs Gamma-ray correlation 

Visual inspection of variability curves in Fig. |2] shows an ap- 
parent correlation between the various flux density peaks of the 
GeV light curve and the optical peaks (Ol to 09, except 06). 
The flaring pattern at y-rays is similar to the QPO-like behavior 
observed at optical frequencies. In addition, the long-term vari- 
ability features are also simultaneous at the two frequencies. To 
compare the flaring behavior of the source at optical and y-ray 
frequencies, we plot the normalized weekly averaged optical and 
y-ray light curves on top of each other (see Fig. [T3]top). A con- 
sistent and simultaneous flaring behavior can be seen between 
JD' = 680 to 1200, however the y-ray variability is less corre- 
lated later. In Fig. [T3] (middle), we show a flux-flux plot of the 
weekly averaged y-ray vs optical V-band data. A clear correla- 
tion among the two can be seen, which is confirmed by a linear 
correlation analysis, yielding rp = 0.36 and 99.996 % confi- 
dence level. The correlation is even stronger in part I, for which 
we find rp = 0.66 and 99.9999 % confidence level. Here, we 
have used the weekly averaged optical flux for the analysis, and 
the uncertainty represents variation of the flux over this period. 

In Fig. [T3] (bottom), we show the cross-correlation analysis 
results of the y-ray and optical data trains. We consider three dif- 
ferent cases to investigate the possible correlation and summa- 
rize the results in Table|7] This analysis reveals that the two-year- 
long GeV and optical data trains are strongly correlated with 
each other with no time lag longer than one week. It is also im- 
portant to note that the strength of correlation is higher before 
the end of the 05/G5 flares than after those flares. 



3.2.5. The orphan X-ray flare 

In order to investigate the origin of the X-ray flare (JD' - 1 120- 
1210, Fig.|2]i, we explore the correlation between X-ray photon 
index and flux. We do not see any systematic change in the X- 
ray photon index (Tx-ray) w.r.t. a change in the flux. The X-ray 
photon index vs flux plot over the flaring period between "5-8" 
(see Fig. |2] for labeling) is shown in Fig. [14] (top) and the esti- 
mated correlation coefficient rp is 0.25 with a confidence level 
of 69 %. Thus, as per correlation statistics, the X-ray photon in- 
dex and flux are not significantly correlated with each other. We 
also notice that the flaring amplitude is similar at soft and hard 
X-rays as shown in Fig. [14] (bottom). The percentage fractional 
variability is 22.5 and 25 in the soft and hard X-ray bands, re- 
spectively. The comparable fractional variability implies that the 
X-ray flare is equally attributed to emission from the soft and the 
hard X-ray bands. 

Although the X-ray light curve of the source is the least sam- 
pled one among all the multi-frequency light curves, we notice 
a flare peaking between "5" and "6" [JD' = 1000 to 1200] (see 
Fig.O. However, due to the gap in the observations it is hard to 
determine the exact peak time of the flare. If we consider that the 
maximum in the X-ray light curve (say X6) is close to the peak 
of the flare, then this epoch coincides with a minimum in the 
optical/GeV flux and it is observed ~ 50 days after the 05/G5 
flares (see Fig.|2]i. 
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Fig. 13. Top: The weekly averaged normalized flux at y-ray and optical V band frequencies plotted vs. time. The flux variations at 
these two frequencies seem to have a one-to-one correlation with each other. Bottom left : y-ray vs optical flux. Bottom right: The 
DCF curve of y-ray vs optical V passband flux using a bin size of 10 days in each case. A: using the complete data as shown in Fig. 
121 B: after removing the data covering the duration of the optical flare 06; C: using the data before flare 06. 



The DCFs of the X-ray light curve with y-ray and radio fre- 
quency light curves do not follow any particular trend as there 
are very few observations available in the X-ray band. A formal 
X-ray vs optical DCF curve (Fig.fTsTi shows a peak at a time lag 
= -(60 + 3) days and another peak at (15 + 3) days. The large 
DCF error bars are due to sparse data sampling of the X-ray light 
curve. In the former case, a negative time lag means that optical 
variations lead the X-ray ones, while in the other case the op- 
posite occurs. An overall inspection of the light curves in Fig.|2] 
reveals that the optical flare (05) is observed ~ 55 days earlier 
than the X-ray flare X6, and 06 appears -12 days later. This 
indicates that the X-ray variability is governed by some other 
eff'ect than the major optical/GeV flares (05/G5), which appear 
strongly correlated. 



3.3. Ha6\o Spectral analysis 

In the following sections, we will study the spectral variability of 
S5 0716H-714 during the different flaring episodes with a focus 
on the good spectral coverage in the radio bands. 



3.3.1 . Modeling the Radio Spectra 

The multi-frequency radio data allow a detailed study of the 
spectral evolution of the two major radio flares, R6 and R8. We 
construct quasi-simultaneous radio spectra using 2.7 to 230 GHz 
data. To perform a spectral analysis of the light curves simulta- 
neous data points are needed. This is achieved by performing a 
linear interpolation between the flux density values from obser- 
vations. A time sampling Af = 5 days is selected for the interpo- 



13 



B. Rani et al.: Radio to gamma-ray variability study of blazar S5 0716+714 



L5 




10 20 30 




0.3-2.0 KeV - 
2.0-10,0 KeV- 



1100 1150 1200 1250 1300 

JD [2454000 +] 

Fig. 14. Top: X-ray photon index vs flux at 0.3 - 10 keV. The 
data points in the box belong to a phase of brightening shown 
in the bottom figure. The X-ray photon index of the source is 
almost constant at 2.25 + 0.25 (shown by a dashed line) over the 
flaring period. Bottom: Soft and hard X-ray light curves of the 
source over the period of high X-ray activity. The flaring activity 
is similar in the two X-ray bands. 
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Fig. 15. The DCF curve of X-ray vs optical V passband flux us- 
ing a bin size of 3 days. 



Table 8. Best-fit spectral parameters for the evolution of radio 
flares using a one-component SSA model 



bin 


Time 


s,„ 


y,„ 


a, 






JD' [JD-2454000] 


[Jy] 


[GHz] 






1 


1096-1101 


0.58±0.09 


95.05±21.78 


0.70±0.26 


-1.15±0.61 


2 


1130-1135 


2.45±0.11 


87.26±7.92 


1.26±0.18 


-0.37±0.13 


3 


1150-1155 


4.64±0.14 


84.92±5.27 


1.15±0.11 


-0.40±0.09 


4 


1173-1178 


7.83±0.34 


80.52±4.59 


1.12±0.11 


-0.62±0.12 


5 


1189-1194 


7.57±0.84 


58.96±8.05 


1.37±0.41 


-0.61±0.30 


6 


1197-1201 


6.42±0.45 


55.78±2.90 


1.06±0.12 


-1.48±0.26 


7 


1204-1209 


5.13±0.24 


60.49±2.22 


0.97±0.08 


-1.56±0.18 


8 


1216-1221 


2.52±0.20 


57.39±6.03 


0.70±0.10 


-1.24±0.33 


9 


1226-1230 


3.23±0.17 


97.03±12.60 


0.39±0.08 


-1.14±0.71 


10 


1238-1242 


3.39±0.20 


106.80±17.90 


0.32±0.05 


-1.30±0.69 


11 


1267-1272 


3.57±0.13 


92.24±7.94 


0.43±0.07 


-0.94±0.43 


12 


1273-1278 


3.21±0.15 


87.92± 14.80 


0.58±0.44 


-0.31±0.12 


13 


1283-1288 


3.42±0.11 


I30.70±32.50 


0.67±0.09 


-0.35±0.17 


14 


1290-1295 


4.32±0.13 


115.30±8.14 


0.69±0.05 


-0.71±0.18 


15 


1298-1303 


6.04±0.20 


74.81±4.44 


1.07±0.10 


-0.47±0.09 


16 


1309-1313 


4.48±0.31 


62.35±9.09 


1.07±0.26 


-0.33±0.17 


17 


1318-1323 


2.77+0.08 


45.66+3.06 


1.56±0.29 


-0.29±0.06 


18 


1340-1345 


1.55±0.06 


40.00±5.22 


1.07±0.27 


-0.18±0.09 



lation. We interpolate the data between the two adjacent observa- 
tions to predict the flux if the data gap is not longer than 5 days; 
however for longer gaps we drop such data points. In Fig. [16] 
we report the spectral evolution of the R6 and R8 radio flares. In 
this figure, the flux densities are averaged values over the 5 days 
binning period and the uncertainties in the fluxes represent their 
variation over this period. 

The observed radio spectrum is thought to result from the 
superposition of emission from the steady state (unperturbed re- 
gion) and the perturbed (shocked) regions of the jet. We con- 
structed the quiescent spectrum using the lowest flux level dur- 
ing the course of our observations. Emission from a steady jet is 
better characterized by a relatively flat spectrum, so we choose 
the steepest one observed on February 17, 2008. The quiescent 
spectrum is shown in Fig.[T6](b). The flux densities were fitted 
by a power law F(v) = Q(v/GHz)"" with Cg = (0.92 + 0.02) Jy 
and ag = -(0.062 + 0.007). 

We fitted the radio spectra using a synchrotron self-absorbed 
spectrum. A syn chrotron self-absorbed (SSA) spectrum can be 
described as (see: lTurler et al.l l2000l iFromm et al.i 1201 1[ for de- 
tails) : 



" l-exp(-T,„(y/v„,)"°-"') 



(2) 



1 - exp (-T„) 

where t„, ^ 3/2 ^ -^1 - §^ - 1 j is the optical depth at the 

turnover frequency, S ,„ is the turnover flux density, v,„ is the 
turnover frequency and a, and q-q are the spectral indices for the 
optically thick and optically thin parts of the spectrum, respec- 
tively (S ~ v"). 

For the spectral analysis, we first subtract the contribution 
of the quiescent spectrum from the data and then used eq. |2]for 
fitting. The uncertainties of the remaining flaring spectrum are 
calculated taking into account the errors of the interpolated data 
points and the uncertainties of the quiescent spectrum. We tried 
two independent approaches to model the radio spectra: (i) a 
one-component SSA model, (ii) a two-component SSA model. 

One-component SSA model : During the fitting process we 
allowed all four parameters {Sm, Vm, at, Q^o) (see eq.|2]i to vary. 
In Fig. [16] (a) we show the 230 GHz light curve with labels 
("numbers"), marking the time of best spectral coverage, for 
which spectra can be calculated. A typical spectrum (for time 
bin "4") is shown in Fig. [16] (c). In Table [8] we list the spectral 
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Fig. 16. The Evolution of the radio spectra: (a) 230 GHz light curve showing different periods over which the spectra are constructed, 
(b) Quiescent radio spectrum; (c) Results of a single component spectral fitting at time bin "4", the dotted line corresponds to the 
quiescent spectrum, the dashed one to the flaring spectrum and the solid line to the total spectrum, (d) The same spectrum fitted 
by a two-component synchrotron self-absorbed model, with the green dashed line showing the individual components and the blue 
solid line showing a combination of the two. A single component model curve is displayed with a dotted-dashed red curve for 
comparison, (e) & (f) The time evolution of S ,„ax vs v,„ax for the R6 and R8 radio flares. The spectral evolution extracted using a 
single-component model is shown by blue symbols and the red symbols denote a two-component model (see text for details). 



parameters of the one-component SSA fit for all spectral epochs 
(bin 1 to 18). In a homogeneous emission region, the spectrum 
is described by characteristic shapes ly oc v^^^ and ly oc v"('"')/2 
for the optically thick and thin domain (s is the power law 
index of the relativistic electrons), respectively. Thus, the 
theoretically expected value of the optically thick spectral index, 
a, is 2.5. While fitting the spectra with a single-component 
SSA model, we find that a, varies between 0.32 to 1.56. This 
deviation of a, from 2.5 indicates that the emission region is 
not homogeneous, and it may be composed of more than one 



homogeneous components. We also notice that the radio spectra 
over the period between the two radio flares R6 and R8 (from 
bin9 to bin 12) can not be described by such a spectral model 
at all. Apparently, these spectra seem to be composed of two 
different components, one peaking near 30 GHz (low-frequency 
component) and the other one at ~ 100 GHz (high-frequency 
component). Consequently, we consider a two-component 
model. 
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Table 9. The best fitted spectral parameters over the evolution of radio flares using a two component SSA model 



bin Time v,„/* 5,,,/ oqI v„,h' S „,h a^h 



JD' [JD-2454Q00] GHz Jy GHz Jy 



1 


1096- 


1101 


20±0 


0.26±0.10 


-0.47+0.59 


98.52+32.74 


0.43+0.13 


-0.51+0.59 


2 


1130- 


1135 


20±0 


0.41±0.16 


-0.12+0.13 


90.00+15.25 


1.98+0.11 


-0.12+0.03 


3 


1150- 


1155 


20±0 


1.37±0.14 


-0.32+0.08 


86.95+6.14 


3.70+0.13 


-0.23+0.05 


4 


1173- 


1178 


20±0 


2.31 ±0.29 


-0.55+0.20 


82.41+4.79 


6.63+0.39 


-0.41+0.08 


5 


1189- 


1194 


20±0 


2.17±1.00 


-0.40+0.47 


59.90+8.55 


6.33+1.05 


-0.62+0.31 


6 


1197- 


1201 


20±0 


2.17±0.24 


-0.39+0.12 


57.01+1.62 


5.57+0.50 


-0.72+0.31 


7 


1204- 


1209 


20±0 


1.46±0.26 


-0.21+0.11 


58.66+2.18 


4.55+0.55 


-0.68+0.38 


8 


1216- 


1221 


20±0 


1.13+0.13 


-0.13+0.05 


57.88+3.02 


2.09+0.60 


-0.91 + 1.50 


9 


1226- 


1230 


18±1 


0.95+0.29 


-1.06+0.69 


128.20+7.33 


3.34+0.10 


-0.46+0.11 


10 


1238- 


1242 


18±1 


1.20+0.41 


-0.76+0.38 


126.10+8.48 


3.40+0.13 


-0.55+0.17 


11 


1267- 


1272 


22±1 


0.76+0.13 


-1.52+0.67 


124.40+4.18 


3.60+0.07 


-0.38+0.03 


12 


1273- 


1278 


23±1 


0.94+0.19 


-1.83+0.98 


129.90+8.47 


3.29+0.12 


-0.33+0.04 


13 


1283- 


1288 


20±0 


1.55+0.18 


-0.16+0.06 


116.30+12.20 


2.24+0.17 


-0.20+0.02 


14 


1290- 


1295 


20±0 


1.65+0.17 


-0.11+0.03 


113.70+15.64 


2.82+0.18 


-0.28+0.13 


15 


1298- 


1303 


20±0 


2.42+0.16 


-0.19+0.03 


75.51+3.52 


4.10+0.17 


-0.41+0.06 


16 


1309- 


1313 


20±0 


2.61+0.41 


-0.47+0.75 


80.93+9.89 


3.27+1.88 


-0.40+0.28 


17 


1318- 


1323 


20±0 


1.82+1.07 


-0.33+0.40 


50.10+5.70 


1.27+1.51 


-0.21+0.27 


18 


1340- 


1345 


20±0 


0.87+0.05 


-0.05+0.01 


55.70+7.80 


0.77+0.05 


-0.26+0.08 



* ; Index I is for low-frequency component and h is for high-frequency component. 



Two-component SSA model : Since the flux densities at cm 
wavelengths are much higher than the extrapolation of the mm- 
flux with a spectral index a, =2.5 for the optically thick branch 
of a homogeneous synchrotron source, we assume that besides 
the mm-submm emitting component, there is an additional spec- 
tral component which is responsible for the cm emission. We 
therefore fit the radio spectra with a two-component model. This 
allows us to fix a,, and set it to 2.5 for both of the compo- 
nents. We also fix the peak frequency of the lower frequency 
component to 20 GHz, as the low-frequency v„, varies between 
1 8 - 25 GHz and the fitting improves only marginally if we al- 
low this parameter to vary. Hence, we study the spectral evo- 
lution of the radio spectra by fixing 0,(1^ - a,{h) - 2.5 and 
v,„i - 20 GHz. Such a scenario appears reasonable and is mo- 
tivated by th e idea of a synchrotron self -absorbed "Blandford- 
Konigl" jet ( iBlandford & Konigl Il979t) and a more variable 
core or shock component. The fitted spectrum using this re- 
stricted two-component model is shown in Fig. [16] (d) and the 
best fit parameters are given in Table|9] A variable low-frequency 
component provides a better fit over bin 9-12. Therefore, we 
consider that both the low- and high- frequency components are 
varying over the time period between the two flares. The two- 
component SSA model describes the radio spectra much bet- 
ter than a single-component model. We therefore conclude that 
the radio spectra are at least composed of two components, one 
peaking at cm wavelengths and the other at mm-submm wave- 
lengths. 

3.3.2. Evolution of radio flares 

In the foll owing we adopt a model of spectral evolution as de- 
scribed bv lMarscher & Gear! (Il985b which considers the evolu- 
tion of a traveling shock wave in a steady state jet. The typical 
evolution of a flare in turnover frequency - turnover flux density 
(S„, - v„,) plane can be obtained by inspecting the R (radius of 
jet)-dependence of t he turnover frequenc y, v,„, and the turnover 
flux density, S ,„ (see lFromm et al.L|201 ll for details). During the 
first stage, Compton losses are dominant and v,„, decreases with 
increasing radius, R, while 5,,,, increases. In the second stage. 



where synchrotron losses are the dominating energy loss mech- 
anism, the turnover frequency continues to decrease while the 
turnover flux density remains constant. Both the turnover fre- 
quency and turnover flux density decrease in the final, adiabatic 
stage. We studied the evolution of the radio flares using the re- 
sults obtained from both one- and two-component SSA models 
and in each case, we obtained similar results. The evolution of 
the R6 and R8 flares in S,,, - v„, plane are shown in Fig. [T6](e) - 
(f). 

In a standard shock in jet model, S,,, oc (iFromm et al.L 

1201 UlMarscher & Geailll985l) where e,- depends upon the varia- 
tion of physical quantities i.e. magnetic fiel d (B), Doppler factor 
(6) and energy of relatiyistic electrons (see iFromm et al.l 1201 ll: 
iMarscher & Gearlll985[ for details). The estimated e, values for 
both the one- and two-component SSA models are given in Table 

m 

IMarscher & Gear (Il985b predi cted a value of ecompton 



-2.5 and IFromm et al 



1: low-frequency component, h: high-frequency component 



. _ (1201 ll) obtain -1.21, whereas 

iBiomsson & AslaksenI ( 2000l) obtained ecompton = -0.43 using a 
modified expression for the shock width. The estimated ecompton 
value for the R8 flare lies between these values, while for the 
R6 flare it is too high to be explained by the simple assump- 
tions of a stan dard shock-in-jet inodel ( see Table [TOl i. For the 
adiabatic stage IMarscher & Gearl (1 19851) derived an ex ponent 
^adiabatic = 0.69 (assuming i = 3) and Fromm et al.l (1201 ll) found 
Cadiabatic = 0.77. We obtain eadiabatic ~ 2 for the R8 flare and ~ 10 
for the R6 flare which is again too steep. The spectral evolution 
of the R8 radio flare can be well interpreted in terms of a standard 
shock-in-jet model based on intrinsic effects. The rapid rise and 
decay of S„, w.r.t. v,„ in the case of the R6 (see Fig.lTSIl flare rule 
out these sim ple assumptions of a stan dard shock-in-jet model 
considered bv lMarscher & GeaH(ll985l ) with a constant Doppler 
factor (6). 

We argue that the spectral evolution of the radio flare, R6 (in 
S m - Vm plane) need to be investigated by considering both the in- 
trinsic variation and the variation in the Doppler factor (i5) of the 
emitting region. lOian et al ] (fT996l) studied the intrinsic evolution 
of the superluminal components in 3C 345 with its beaming fac- 
tor variation being taken into account with a typical variation of 
the viewing angle by 2 - 8°. In the s tudy of the spectr al evo- 
lution of the IR-mm flare in 3C 273, lOian et all (1 19991) found 
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that the bulk acceleration of the flaring component improves the 
fit of the spectral evolution at lower frequencies. Therefore, it 
is worthwhile to include a variation of 6 along the jet axis in 
our model, which we parametrize as 5 oc Such an approach 
could easily explain the large variation in the observed turnover 
flux density, while the observed turnover frequency kept a ne arly 
constant value or changed only slightly (iFromm et al.Ll20ll1) . 

We consider the evolution of radio flares in the framework 
of dependencies of physi cal parameters a, s and d following 
iLobanov & Zensusid 19991) . Here, a, s and d parametrize the vari- 
ations of (5 oc R'^, B oc R " and N{y) oc y^ ' along jet radius. Since 
it is eviden t that e values do not difi'er much for different choices 
of a and s (ILobanov & Zensus, 1999) , we assume for simplicity 
that s w constant and for two extreme values of a = 1 and 2, we 
investigate the variations in b. The calculated values of b for dif- 
ferent stages of evolution of radio flares are given in Table [TOl It 
is important to note that the Doppler factor varies significantly 
along jet radius during the evolution of the two radio flares. 

Moreover, the turnover frequency between the Compton and 
synchrotron stages (v, ) and the synchrotron and adiabatic stages 
(vy) in the S„, - v,„ plane characterize th e observed behavior 
of the radio outbursts (IValtaoia et al 1 119921) . In a shock induced 
flare, the shock strength reaches its maximal development at v, 
and the decay stage starts at Vf. In Fig. [16] (e) - (f), we display 
by dashed lines the frequencies Vr and vj. The shock reaches its 
maximal development at 80 GHz for the R6 flare and at 74 GHz 
for the R8 flare. The observed behavior of the outburst depends 
on V, . In a shock induced flare, the observed frequency {v„bs) is 
less than v, in the case of low-peaking flar es, while Vobs > Vr 
for high-peaking flares (IValtaoia et al 1 119921) . We therefore con- 
clude that both the R6 and R8 radio outbursts are low-peaking 
radio flares and are in quantitative agreement with the formation 
of a shock and its evolution. 

3.3.3. Synchrotron Spectral Break 

The source was observed at IR frequencies with 
the Spitzer Space Telescope on December 06, 2007. 
We obtained the IRAC+MIPS photometric mea- 
surements at 5 - 40 pm from the Spitzer archive 
(http://sha.ipac.caltech.edu/applications/Spitzer/SHA/). Since 
the source has been observed at radio wavelengths over this pe- 
riod, we combine the cm - mm and IR observations to construct 
a more complete broadband synchrotron spectrum. The com- 
bined radio - IR spectrum is shown in Fig. [17] The red curve rep- 
resents the best fitted synchrotron self-absorbed spectrum with a 
break at a frequency of vj, = (1.3 + 0.1) x 10"* GHz. The best-fit 
parameters are: S,„ = (1.03 + 0.02) Jy, v„, = (45.74 + 3.12) GHz, 
at - (0.33 + 0.01) and the spectral index of the optically 
thin part (cq) is -(0.38 + 0.09) and -(0.66 + 0.07) above and 
below the break, respectively. Hence, modeling of the radio 
- IR spectrum provides strong evidence for a break in the 
synchrotron spectrum at v;, ~ 1.3 x lO'* GHz with a spectral 
break Aa - 0.28 + 0.1. We can also include the optical V 
passband flux from the AASVO (see section 12.2b to estimate 
the spectral break. This leads to a steeper spectral index with 
ao = -0.88 + 0.03 and a break Aa = 0.51 + 0.09. 

The spectral break could be attributed to synchrotron loss 
of the high energy electrons. It is widely accepted that syn- 
chrotron losses result in a steepening of the particle spectrum 
by one power and a steepening of the emitted synchrotron spec- 
trum by a half-power (Reynolds, 2009; Kardashey, L962)- Also, 
synchrotron-loss spectral breaks differing fro m 0.5 could be pro - 
duced naturally in an inhomogeneous source (iRevnoldsl l2009l) . 
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Fig. 17. Radio-IR spectra using Spitzer observations. The red 
curve is the best fitted synchrotron self-absorbed spectra with 
abreakat(1.3+0.1)xlO''GHz with a spectral break Aa - 0.28 + 
0. 1 . The green line represents the spectral fitting including opti- 
cal data point and this leads to spectral break Aa = 0.51 + 0.09. 

As Vb is mainly determined by synchrotron loss, it depends on 
the magnetic field strength. One can estimate the minimum- 
energy magnetic field strength using the following relation given 
byllJeavens & Meisenheimer ( 1983): 

Bbreak = 2.5 x IQ-'/i]^' L-^'\-^'^' G (3) 

where ySi .c is the speed of the upstream gas related to the shock, 
L is the length of the emission region in kpc (at v < v;,) and v;, 
is the break frequency in GHz. For relativistic shocks /3i is close 
to 1 . We constrain the length of the emission region L using the 
variability timescales at 230 GHz as this is the closest radio fre- 
quency to Vb- Using L < 0.04 x 10""^ kpc (see Section l3.4.3b . we 
found Bbrmk ^ 0.14 G. The minimum energy condition implies 
equipartition of energy, which means Bbreak ~ B^g (equipartition 
magnetic field). 

3.4. Brightness Temperature, size of emission region and Jet 
Doppler factors 

3.4.1. Brightness Temperature T"/'' 

The observed rapid variability implies a very compact emission 
region and hence a high brightness temperature if the variations 
are intrinsic to the source. Assuming a spherical brightness dis- 
tribution for the variable source and that the triggered flux varia- 
tions propagate isotropically through the source, then the light 
travel time argument implies a radius d < cAt for the emis- 
sion region where At is the time interval of expansion. So, the 
flux variability observed in radio bands allows us to estimate 
the brightness temperature of the source usi ng the relation (see 
lOstorero et all 120061: iFuhrmann et"an l2008l for details). 

Tt"" = 3.47 X 10^ ■ AS, I ^^^f K (4) 
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Table 10. Different states of spectral evolution and their characteristics 



Flai"e 


Time 


bin 




6 


b 


Stage 




JD' [JD-2454000] 




(1 component SSA) 


(2 component SSA) 


s=2.2, a=l-2 




R6 


1096-1178 


1-4 


-7±3 


-8±3 


0.7 


Compton 




1178-1194 


4-5 








-0.07 


Synchrotron 




1194-1221 


5-8 


I0±2 


11±3 


2.6 


Adiabatic 


R8 


1283-1303 


13-15 


-0.9±0.1 


-1.2±0.2 


0.4 


Compton 




1298-1345 


15-18 


1.8±0.2 


2.5±0.5 


-2 


Adiabatic 



d cc r'' , B oc R-" and N(y) cc y-' 



Table 11. Variability Brightness temperatures 
R6 flaie 



Frq. (GHz) 


t" 

var 






e 


(GHz) 


(days) 


(10^' K) 




(mas) 


15 


61 


154 


10 


0.091 


37 


65 


62 


7 


0.068 


86 


60 


13 


3 


0.027 


230 


50 


3 


2 


0.015 


R8 Flare 


15 


37 


307 


14 


0.077 


37 


18 


109 


9 


0.025 


86 


25 


55 


7 


0.021 


230 


10 


5 


3 


0.004 



a: the lower value corresponds to the rising phase of flare while the 
higher to the decaying phase. 



and field energy dScott & Readheadi [19771) : Te.eg ~ 5 x 10"' K 
which is derived on the basis of an argument that this limit bet- 
ter reflects the stationary state of a synchrotron sour ce which for 
many sources yields Tb < 10" K ( iReadheadll 1994 e.g.). In this 
case, the calculated Doppler factor values using the equiparti- 

tion limit, 6yar,eq - (1 +z) ^*^Tg'''' /5 x 10'°, become higher by a 

factor of 4.47 i.e. 6yar,eq - 4.47 X d^.ar- 



3.4.3. Size of the emission region 6 

One can obtain the size of the emission region using the calcu- 
lated Doppler factors (6) and variability time scales {t„r) '■ 



0.173^5(1 +z) mas 



(6) 



where A5 ^ is the change in flux density ( Jy) over time tvar,A 
(years), di is the luminosity distance in Mpc, A is wavelength in 
cm and z is the redshift of the source. Here and in the following 
calculations we will use z - 0.31, which yields a l uminosity dis- 
tance, di - 1510 Mpc (see lFuhrmann et al.[|2008l for details). 

Two major outbursts (R6 and R8) are observed in the source 
at 15 GHz and at higher radio frequencies. To calculate the 
brightness temperature, we have used the rising time of the flares 
(see Table |5]) separately for the two flares at 15, 37, 86 and 
230 GHz, as these are the best sampled light curves. The radio 
flares follow a slow rising and fast decaying trend, so we calcu- 
late t„,r for both the rising and the decay phases of the flares. The 
calculated for the two radio flares are listed in column 2 of 
Table nn and the apparent brightness temperatures (1'°'^'') are in 
column 3. 



3.4.2. Doppler factor from variability timescales 6,.^, 



The calculated is one to two orders of magni- 

tude higher than the IC- li mit T 'J;"/;, of Tb ~ 10'^ K 
jKellermann & Paulinv-TotR Il969h at all frequencies up to 
230 GHz. We notice that Tg'''' exceeds Tg""^ even at short-mm 
bands. The excessive brightness temperature can be interpreted 
by relativistic boosting of the radiation, which gives to a lower 
limit of the Doppler factor of the emitting region 



(l+z) 



r"PP 



101 



(5) 



Here a is the spectral index of the optically thin part of the radio 
spectrum. We obtained a,i,i„ — -0.23 to -0.91 for the R6 radio 
flare and a,/,,,, = -0.20 to -0.41 for the R8 flare (see Section 4.1 
for details). The calculated 5,,,,^ values are listed in column 4 of 
Table fm We obtain dn,,- > 14 for the two radio flares. 

In addition, we can also use the intrinsic brightness temper- 
ature limit based on the equipartition between particle energy 



The angular size calculated using dmrjc are listed column 
5 of Table [TT] We obtain that the estimated value of the angular 
dimension lies between 0.004 - 0.09 mas. Again, 6 will be a 
factor of 4.47 higher if we use 6yar,eq- In linear dimensions, the 
size of the emission region ranges between (0.6 - 12.1)xl0'^ 
cm. 



3.4.4. Inverse Compton Doppler factor 5ic 

One can constrain the Inverse Compton Doppler factor (5/c) by 
comparing the expe cted an d observed fluxes at high energies 
(see Ghiselli ni et all I1993L for details). The IC Doppler factor 
is defined as 



5,c^{mS,„{\+z)\ 



5y0«'-4,v)^(5-3a) 



1/(10-60-) 



(7) 



where Vc is the synchrotron high frequency cut-off in GHz, Sm 
the flux density in Jy at the synchrotron turnover frequency v„,, 
S jc the observed y-ray flux in Jy (assumed to arise from the 
IC process) at Vy in keV, a is the spectral index of the opti- 
cally thin part ol^ the spectrum, 0y the source's angular size in 
mas and f{a) ^ 0.14 - 0.08q'. The apparent variability size is 
calculated using eg. [6] For the high energy cut-off we follow 
iFuhrmann etalTdlOOSi) and use ~ 5.5 x 10^ GHz. 

For these calculations, we used equals 9 days, which 
is the fastest variability timescale for the R6 flare at 86 GHz. 
a is obtained from the SSA modeling (see Section 13.3. Il l and 
the estimated values are given in Table fT2] die is calculated for 
the same four time bins which we use to model the broadband 
SEDs of the source (see Section l33T l. The estimated values for 
die are given in Table fT2l We find that duiing the four different 
activity states of the source Sjc > 20. 



18 



B. Rani et al.: Radio to gamma-ray variability study of blazar S5 0716+714 



Table 12. Brightness temperature 



Time bin 6ic Used parameters 



Binl 


SiC.O.SKeV 


> 11 


SiQ = 3.71 Jy, a = 


-0.74, 








So.sKev =2.97x10- 


Jy 




SlCJ.SKeV 


> 20 


5go = 3.71 Jy, a = 


-0.74, 








Slo.SKeV =8.27x10 


Jy 




SlC.lOQMeV 


> 14 


5go = 3.71 Jy, a = 


-0.74, 








S looMev =1.17x10 


Jy 


Bin2 


SiC.Q.iKeV 


> 11 


540 = 1.68 Jy, a = 


-0.52, 








So.sKev =2.97x10- 


*^ Jy 




SiC.l.SKeV 


> 14 


540 = 1-68 Jy, a = 


-0.52, 








SrsKev =4.57x10- 


Jy 




SiClOOMeV 


> 14 


540 = 1-68 Jy, a = 


-0.52, 








5 looMev =1.50x10 


Jy 


Bin3 


SiC.a.SKeV 


> 14 


582 = 9.89 Jy,Q- = 


-0.76, 








So.sKev =3.85x10- 


* Jy 




SiC.l.SKeV 


> 15 


582 = 9.89 Jy, a = 


-0.76, 








SjsKev =2.97x10" 


Jy 




SlC.lOOMeV 


> 17 


582 = 9.89 Jy, or = 


-0.76, 








SmMev =2.94x10 


Jy 


Bin4 


SiC.O.SKeV 


> 12 


5,8 = 3.85 Jy, a = 


-0.78, 








So.sKev =2.97x10- 


•5 Jy 




SiC.l.SKeV 


> 12 


5,8 = 3.85 Jy, a = 


-0.78, 








Si.sKev =2.97x10- 


« Jy 




SiC.lOQMeV 


> 12 


578 = 3.85 Jy,Q- = 


-0.78, 








5 \ooMev =2.08x10 


Jy 



3.4.5. Gamma-ray Doppler factor 6y 

It is also possible to obtain a limit on the Doppler factor 5 by 
considering that the high-energy y-ray photons can collide with 
the softer radiation to produce e=^ pairs with the assumption that 
the bulk of the high-energy emission (y-rays and X-rays) is pro- 
duced in the same emission region. The cross- s ection of this 
process is maximized at ~ ctjIS (see ISvensso for de- 

tails), where a-j is the Thomson scattering cross-section. This 
leads to a lower limit on S with the req uirement that Tyy{v) < 1 
(iDondi & GhiselliniL[T99llFinke et al.[l2008l) : 

6y > [ ^ '- "° (8) 

where a is the power law index of the synchrotron spectrum i.e. 
ff" oc cTj is the scattering Thomson cross-section, is the 
electron mass, e\ — Ej{meC-') is the dimensionless energy of a 
y-ray photon with energy E for which the optical depth of the 
emitting region Tyy - 1 . For the highest energy GeV (207 GeV) 
photon observed in the source (Rani et al. 2012), we obtain e - 
207GeV/(5.11 x lO-'^GeV) = 4 x 10"* and e"' = 2.4 x 10"^ 
Using /"," = 3.88 x lO"" ergs cm-^ s"', we obtain 5^ > 9.1 
The d etection of the source at above 400 GeV (Anderhub et al.L 
120091) constrains (5^ > 9.8 

3.4.6. Magnetic field from synclirotron self-absorption 

It is also possible to constrain the magnetic field using the 
standard s ynchr otron self-absorption expressions. Following 
iMarsched (Il987l) . an expression for the magnetic field B in a 
homogeneous synchrotron self-absorbed region is given by: 

BsA [G] = lO-'b(a)0'^^„,S-r (y^) , (9) 

where b{a) de pends on t he op tically thin spectral index o-,/,,,, 
(see Table 1 in iMarscheR 1 1 9871) . S,n is the flux density, is the 



source's angular size at the synchrotron turnover frequency Vm 
and 6 is the Doppler factor. The size of the emitting region re- 
sponsible for the observed variations can be constrained using 
m m-VLBI me a surem ents of the core region of S5 0716-1-714 
by iBach et alj (l2006h : < 0.04 mas. Using b{a) = 3.13, 
S„t - 3.89 Jy, Vm = 80 GHz, we calculate a lower limit of 
the magnetic field Bsa in the range of (0.0078-0.0198)5 mG. 
Using 6 > 1 at v„i ~ 80 GHz (see Table [TTI). we obtained Bsa ^ 
0.05 to 0.14 G. The size of the emission region constrained us- 
ing the causality arguments, ~ 0.0027 mas at v,„ ~ 80 GHz 
(see Table fTTTi gives Bsa ^ 0.03 G. These calculations constrain 
Bsa > 0.03- 0.14 G 

3.4.7. Equipartition IVIagnetic field and Doppler factor 

The equipartition magnetic field B^^j, which minimizes the to- 
tal energy Ej^t — (1 -f- k)E^ + Eg (with relativistic particle en- 
ergy Ef, ~ Z?-'^ and energy of the magneti c field Eb ~ B^) , 
is given by the follow ing expression (e.g. iBach et all 120051; 
iFuhrmann et"ani2008h : 

B,, = {4.5- {I +k)f{a,v„Vb)LR'f" (10) 

here k is the energy ratio between electrons and heavy parti- 
cles, L is the synchrotron luminosity of the source given by 
L - 4-Kdj^(l + z) j^" S dv, R is the size of the component in 
cm, S „, is the synchrotron peak flux in Jy, Vm is the synchrotron 
peak frequency in GHz and /(o-, v„, v^) is a tabulated function 
depending on the upper and lower synchrotron frequency 

cutoffs Va,Vh- Using = 10^ Hz, vt - 5.5 • 10'"^ Hz, and 
/(-0.5, 10^ 10") = 1.6 X \Q\ we obtain 

B,, = 531x\Q'^[kS„,v,ndlR-^f' G (11) 

Using Beq > 0.14 (see Section |333]l, S ,„ = 3.89 Jy, v„, = 
80 GHz, R = 2.90 x lO"" - 1.2 x 10"* cm (estimated using 
tvar - 25 days at v„, - 86 GHz), the above expression yields 
k - \. A small value of k implies that the jet is mainly composed 
of electron-positron plasma. 

Equation |9] and [TT] give different dependencies of the 
magnetic field on S, i.e. Bsa ~ (> and B^g ~ ^2/7tr+i rpjjj^ yields 

Beq/ Bsa = ^Iq"- Adopting the above numbers, we obtain 
Doppler factors 6eqB in the range of 14 - 20 (for a - -(0.35 to 
0.7)). 



3.4.8. Comparison of the estimated parameters 

The apparent brightness temperature Tb obtained from the day- 
to-day variations exceeds the theoretical limits by several orders 
of magnitudes. Although Tb decreases towards the mm-bands, it 
is still higher than the IC-Umit (lO'^ K). Tb exceeds 10"^ K at 
15 GHz and 10'^ K at 230 GHz. We have obtained lower limits 
to the Doppler factor of the source using different methods as 
discussed in the earlier sections. These methods reveal a range 
of consistent lower limits to the Doppler factor with S^a,- > 14, 
6ic > 20, 6y > 10 and d^q^B ^ 20. Comparing the Doppler factor 
estimates obtained with different methods seems to suggest that 
6 > 20. An independent approach to estimate 6 is spectral mod- 
eling of the broadband SEDs, and this gives 6 = 25 (see Section 
13.5b . which is in agreement with the former values. These lim- 
its are in good agreement with the estima tes based on the re cent 
kinematical VLBI studies of the source (IBach et all I2OO6I) and 
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the IC Doppler factor limits obtained bv lFuhrmann et all (1200 8h . 
As 5eq,B agrees fairly well with the 6 values derived from the 
other methods, we conclude that the emission region is in a state 
of equilibrium. 

The estimated magnetic field value from the broadband spec- 
tral modeling lies between 0.05 and 1 G. A break in the optically 
thin power-law slope at a wavelength of ~ 23 /im constrains 
the equipartition magnetic field to B^q > 0.36 G. We obtained 
BsA ^ 0.14 G from the synchrotron self-absorption calcula- 
tion. The size of the emission region {ff) derived on the basis 
of causality arguments lies between 0.004 - 0.091 mas which 
agrees fairly well with the size of emission regio n constrained 
using mm-VLBI measurements (iBach etal.Ll2006l) . 

3.5. The Complete spectral energy distribution 

The broadband monitoring of the source over several decades of 
frequencies allows us to construct multiple quasi-simultaneous 
SEDs. The SEDs of the source constructed over 4 different 
periods of time are shown in Fig. [18] These time bins 
reflect different brightness states of the source and each time 
bin has a width of 10 days. The variation in flux over the bin 
width is shown as error bars in the SED plots. We construct the 
broadband SEDs for the following activity periods ; 
BinlB : Radio-mm(steady), optical(high), X-ray(steady), 
GeV(low) 

Bin2 : Radio-mm(low), optical(flaring), X-ray(low), 
GeV(flaring) 

Bin3 : Radio-mm(flaring), optical(flaring), X-ray(flaring), 
GeV(low) 

Bin4 : Radio-mm(steady), optical(low), X-ray(steady), 
GeV(steady) 

The double-humped structures of the broadband SEDs can 
usually be modeled by both leptonic and hadronic models (e.g. 
iBoettcher et al.L l2012h . Here we have used a quasi-equilibrium 
version of a leptonic one-zone jet model as described by 
IBoettcher et al.l (l2012h . In this model, the observed radiation is 
assumed to be originating from the ultra-relativistic electrons 
(and positrons) in a spherical emission region of co-moving ra- 
dius Rb propagating with relativistic speed fSrc (F is bulk Lorentz 
factor) along the jet, which is offset by an angle 6 w.rt the line- 
of-sight. We fix 6 to be such that the bulk Lorentz factor, F 
equals the Doppler factor, 6, which, for highly relativistic motion 
(F » 1) implies = 1/F. The emitting electrons are assumed to 
be instantaneously accelerated into a power-law distribution of 
electron energy, E^, - ynieC^, of the form Q{y) - 2oy with q 
being the injection electron spectral index and yi and 72 are the 
low- and the high-energy cutoffs. 

An equilibrium in the emission region between particle in- 
jection, radiative cooling, and escape of particles from the emis- 
sion region yields a temporary quasi-equilibrium state described 
by a broken power law. The particle escape is parametrized 
through an escape timescale parameter ^ > 1 so that tgsc - riR/c. 
The balance between the particle escape and radiative cool- 
ing will lead to a break in the equilibrium particle distribu- 
tion at a break Lorentz factor y/,, where fg^c = tcnoAy)- The 
cooling timescale is calculated self-consistently taking into ac- 
count synchrotron, SSC and EC cooling. Depending on whether 
jb is greater than or less than y\, the system will be in the 
slow cooling or fast cooling regime, respectively, leading to dif- 



we use 'bin' for radio spectra and 'Bin' for radio to GeV spectra 



ferent spectral indices of the equilibrium electron distribution 
iBottcher & Chiand (l2002l) . 

In the fitted model the number density of injected particles is 
normalized to the resulting power in ultra-relativistic electrons 
propagating along the jet given by, 

yn{j)dy (12) 

The magnetic field is considered as a free parameter in 
the emission region. The Poynting flux along jet is Lb - 
nRlr~fircuB, where ub - B- l(%n) is the magnetic field energy 
density. The equipartition parameter eg - Lb/L,, is calculated 
for each fitted model. 

After evaluation of the quasi-equilibrium particle distribu- 
tion in the emission region, our code calculates the radiative 
output from the synchrotron, SSC, and EC emissions self- 
consistently with the radiative cooling rates. The external radia- 
tion field, which serves as seed photons for EC scattering, is as- 
sumed to be isotropic in the stationary AGN rest frame. Its spec- 
trum can be chosen to be a thermal blackbody with temperature 
V, and radiation energy density Ue„, or a line-dominated spec- 
trum (or a combination of the two). The direct emission from 
this external radiation field is added to the emission from the jet 
to yield the total modeled SED, which we fit to the observations. 

We first tried to fit the SEDs with a pure SSC model, as this 
has fewer free parameters than the EC version of the model. 
However, except for the SED of bin 1 (see below), pure SSC 
models typically fail to reproduce the Ferm//LAT spectra of the 
SEDs. Also a detailed study of the y-ray spectrum of the source 
(Rani et al. 2013), shows a significant correlation between de- 
tection of the high energy GeV photons and change in spectral 
slope below and above the break energy, which suggests that 
BLR opacity has a significant impact on the observed spectral 
breaks. Therefore, we included an external radiation component, 
as outlined above, to produce SSC+EC fits. 

The fitted models are shown in Fig. [TS] and the best-fit pa- 
rameters are given in Table [13] The pure SSC model does a 
moderately good job in describing the SEDs of the low states, 
although the y-ray spectra appear systematically too steep. The 
SED of Binl is well fitted with the SSC model, while for the 
other time bins an EC component is required to fit the GeV spec- 
tra. The high-state is very problematic for the SSC model as it 
would require a much lower magnetic field (far below equiparti- 
tion) and - in the case of Bin 2 - a very large emission region, in 
conflict with the often observed intraday optical variability. All 
the low-state fits are possible with parameters close to equipar- 
tition between relativistic electrons and the magnetic field. A 
model including external Compton generally does a better job 
in reproducing the entire SEDs (including the y-ray spectrum), 
if one uses an external radiation field dominated by Ly-a emis- 
sion from a putative broad line region (BLR). For S5 0716+714, 
we found that the radiation field energy density of this external 
field varies between 10"* to 10"^ ergs cm"-', which is a factor of 
~1000 lower than what we expect for a typical quasar. However, 
this is a reasonable value for a BL Lac like S5 0716+714 which 
is known to have a featureless spectrum. Furthermore, this low 
BLR energy density value naturally explains the origin of y-ray 
spectral breaks observed in the source. Moreover, the low BLR 
energy density is consistent with the non-detection of emission 
lines. Parameters close to equipartition can be used for all time 
bins, including the high states. 

At first glance the fits look good, but in more detail the fit to 
the radio data for some bins is relatively poor In the EC model, 
the model fits the cm-radio data quite well, but is much below 
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Table 13. Parameters of SSC and EC fits to SED of S5 0716+714 



Parameters 


ami [JD =«45-s55J 


Bin2 [JD' = 


1 1 1 n 11 ^ on 

1 1 10-1 120J 


Bin3 [JD' = 


lion 11 rvm 

1 180-1 190J 


Bin4 [JD' 


1 1 n 1 m 

= 1210-1220] 




SSC 


SSC 


EC 


SSC 


EC 


SSC 


EC 


n 


l.SxlO"* 


LlxlO"* 


4.0x10'' 


2.5x10"* 


1.8x10"* 


3.0x10"* 


2.5x10"* 




l.OxlO'' 


2.6x10"^ 


6.5x10"^ 


2.0x10^ 


2.0x10' 


l.OxlO'' 


l.OxlO' 


q 


3.10 


3.20 


3.40 


3.15 


3.10 


3.45 


3.45 


1 


25 


100 


25 


25 


25 


25 


25 




1 

i 


U.UJ 


n 7 


n Q 


yj.yj 


U.o 


1 

1 


r 


25 


25 


25 


25 


25 


25 


25 


Ri, (cm) 


1.25x10"^ 


1.7x10" 


2.0x10'* 


1.4x10'* 


2.0x10'* 


7.5x10" 


7.5x10" 


9 (degree) 


2.29 


2.29 


2.29 


2.29 


2.29 


2.29 


2.29 


LAW^] (ergs-') 


1.33 


26.99 


4.15 


4.298 


4.31 


3.09 


2.48 


es 


1.61 


0.063 


1.11 


0.87 


1.59 


0.27 


0.53 


Tf vi K 






Ly-a 




Ly-a 




Ly-a 


E„., (ergcm-^) 






1.7x10"' 




3.0x10"* 




1.0x10"' 



7i , 72 : Low- and High-energy cutoff 

q : Injection electron spectral index 

T] : Electron escape timescale parameter 

B (G) : Magnetic field at z=0 

F : Bulk Lorentz factor 

Ri, (cm) : Blob radius 

6 (degree) : Observing angle 
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Fig. 18. Broad band SEDs of S5 0716+714. Each SED is constructed using 10-day averaged multi-frequency data. The error bars 
represent the variation of flux over 10 days in each bin. Pure SSC models are shown as thick dashed curves. For EC fits, the total 
model SEDs are the thick solid curves; the synchrotron components are dotted, the SSC components are dot-dashed, and the EC 
components are thin dashed curves. 
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the mm data for Bin3. The model for Bin4 does not fit the radio 
data at all (see Fig. [18]). So, in general the model under-predicts 
the radio flux at mm and cm bands. This indicates the possibil- 
ity of a missing spectral component at cm-mm wavelengths. We 
have seen in Section 3.3.1 that a two-component model better 
describes the radio spectra. Therefore, we conclude that an ad- 
ditional synchrotron component is required to fit the broadband 
SED at mm to cm wavelengths. 



4. Discussion 

The densely sampled multi-frequency observations of the BL 
Lac object S5 0716-1-714 over the past three years allow us to 
study its broadband flaring behavior and probe into the physi- 
cal process, location and size of the emission regions. We found 
a direct connection between GeV and optical flares, and major 
flares propagate down to radio wavelengths. The radio outbursts 
seem to be smeared out at 10 GHz and lower frequencies. An or- 
phan X-ray flare lags the major optical-GeV flare (05-G5) by 
~ 55 days and the X-ray emission is produced by both syn- 
chrotron and inverse-Compton mechanisms. It seems that the in- 
teraction of shocks with the underlying jet structure might be 
responsible for optical and high energy emission, and opacity 
plays a key role in the time-delayed emission at radio wave- 
lengths. 



4.1. Broadband correlation alignment 

Following the broadband cross-correlation analysis (13.2) . we 
plot the estimated time lag as a function of frequency in Fig.[T9l 
Fig. [19] (a) shows the plot of the time lag measurements at dif- 
ferent frequencies for the R6 flare using 15 GHz as the reference 
frequency (see Fig. [2]i. It has become evident in Section 13.2. II 
that the time lags (w.rt. 15 GHz) increase with frequency and 
follow a power-law as a function of frequency with a slope ~ 0.3. 
If we extend the fitted power law to optical frequencies, then the 
R6 flare meets the 05 flare, which is observed ~ 60 days earlier 
than the R6 flare. A formal cross-correlation between optical and 
radio frequency light curves indicates a significant correlation 
with a delay of ~ 60 days at radio wavelengths. The dashed line 
in Fig. [19] (a) connects the simultaneous optical - GeV flares. 
The optical-GeV correlation shows no time lag among the flares 
at the two frequencies, i.e. 05 correlates with G5 and 04 with 
G4; but there is no respective GeV flare for 06. The nearby op- 
tical. X-ray and GeV flares are shown with their possible time 
lags w.r.t. R6. The allowed time range of the peak of the X-ray 
flare is marked with an arrow. In Figure[T9](b), we show a similar 
plot for the R8 flare. Both of these figures provide a one-to-one 
connection of the broadband flares based on our analysis. 



(a) 
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Fig. 19. (a) A plot of time lag measurements vs frequency using 
15 GHz as the reference frequency for the radio flare R6. The 
best fitted power law at radio and mm frequencies is extended 
up to the optical wavelengths. The near by optical. X-ray and 
GeV flares are shown with their possible time lags w.rt. R6. (b) 
A similar plot for the R8 flare (see Fig. [2] for flare labeling) In 
both plots, the dashed lines indicate the SSC process with simul- 
taneous optical-y-ray events. 



4.2. Origin of Optical variability 

During our observations, the source was highly active at opti- 
cal frequencies showing multiple flares roughly separated from 
each other by ~ 60 — 70 days, superimposed on a long-term 
variability trend at a ~ 350 day timescale. The periodogram 
analysis reveals two significant peaks at ~ 63 and 359 day 
timescales. A more robust analysis using the power spectrum 
density method implies that the significance of a detection of 
a quasi-periodic signal at the frequency corresponding to these 
timescales is only 2 cr. Thus, the significance of detection re- 
mains marginal. However, it is important to note that periodic 



variations at a year timescale has also been observed earlier in 
flie source (iRaiteri et al.ll2003l) . 

During the two years of our observations, we found that the 
long-term variability amplitude of the source remains almost 
constant at about 1.3 magnitudes. A constant variability ampli- 
tude can be int erpreted in terms o f variations of the Doppler 
boosting factor (IRaiteri et al.Ll2003l) . The change in 6 can be due 
to either a viewing angle (6) variation or a change of the bulk 
Lorentz factor (T) or maybe a combination of both. We notice 
that a change in (5 by a factor of 1 .2 can be easily interpreted as a 
few degree variation in 6, while it requires a noticeable change of 
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the bulk Lorentz factor. We therefore propose that the geometry 
significantly affects the long-term flux base-level modulations. 
Such variations are very likely orig inating as a relat ivistic shock 
traces a spiral path through the jet (iMarscherl Il996l) . 



4.3. Origin of y -rays 

The source displays substantial activity at y-rays during the high 
optical activity period. This is to be expected in leptonic mod- 
els, as the same electrons radiating the optical synchrotron pho- 
tons would emit y-rays through the inverse Compton scatter- 
ing process. Here, we observe a similar flaring behavior at the 
two frequencies. We also found that the flux variations at opti- 
cal and GeV frequencies are significantly correlated with each 
other (on weekly timescales) and corresponding to each optical 
flare "Ol" to "09" (except 06) there is a local maximum "Gl" 
to "G9" at GeV frequencies. In addition to that the variability 
timescales (both short and long) are also comparable for optical 
and GeV light curves. We note that the ratio between the high 
and low y-ray flux levels is about 15, while in the optical band 
the same ratio is of the order of 3.7. Thus, the y-ray flux density 
appears to vary as the square of change in the optical flux density. 
This reflects a quadratic dependence of the GeV flux variations 
compared to optical variability. This favors a SSC interpretation. 
However, we would also like point out that a weak EC contribu- 
tion is also required in order to model the GeV spectrum of the 
source. 



4.4. Opacity and delay at radio wavelengtlis 

The source reaches a maximum in simultaneous optical - GeV 
flaring activity at "5" (see Fig.|2] flares: 05-G5). This maximum 
coincides with the beginning of a major radio outburst "R6" at 
230 GHz. The R6 radio flare is observed ~ 65 days later than 
the optical flare 05 at 230 GHz. The R6 flare is followed by 
another radio flare, R8, with a moderate level of flux activity 
between the two. The two major 230 GHz radio outbursts (R6 
and R8) are smoothed and delayed at lower radio frequencies 
till 15 GHz. The flaring activity seems to be completely washed 
out at < 10 GHz. The estimated time lag (using 15 GHz as ref- 
erence frequency) at each frequency as a function of frequency 
follows a power law with a slope ~ 0.3. Delayed emission at 
lower frequencies is a clear indication of opacity effec ts due to 
synchrotron self-absorption (iKudrvavtseva et aTTl201 ll) . 

As per the cross-correlation analysis, the optical - radio 
variability is found to be significantly correlated, with the 
flux variations at optical frequencies leading those at radio 
bands by ~ 60 days (see Section 13.2.2b . Most earlier studies 
on the radio - optical correlation have shown that the radio 
even ts la g behind the optic al ones by sever al weeks or months 



Tomi koski et al, 1994; IClements et al.LfT99 5: Villa ta et al 



2007 



2011 



iRaiteri et al.L 120031: iJorstad et al.L 120 lOt lAgudo et al 
and references therein). Similar variability timescales (~ 
90 and 180 days) at optical and radio frequencies again hint at 
a co-spatial origin of variability. It is worth pointing out that the 
long term variability timescales are common at optical (and also 
at y-rays) and radio frequencies. The fast repetitive optical/y- 
ray flares are not observed at radio wavelengths. Therefore, it is 
not unreasonable to suggest that the long term variability fea- 
tures observed at optical/GeV frequencies propagate down to 
radio frequencies with a time lag of ~ 60 days. As we notice 
in Section [3. 3.1! the two radio outbursts are low-peaking flares. 
Thus, a 60 day time delay between the optical and radio activ- 



ity em phasizes the optical fl ares being the precursor of the radio 
flares Jvaltaoiaet aL[ll992h . 



4.5. Origin of tine orplian X-ray flare 

When the source is flaring at optical - GeV frequencies, it is 
quiet at X-rays. Although it is hard to locate the exact peak time 
of the X-ray flare, it is obvious that the maximum of the X-ray 
flux peaks ~ 50 days later than the major optical - GeV flares 
(05-G5) (see Fig.|2]i. At this epoch, the source was in a relatively 
steady state at optical/GeV frequencies while there is another 
bright optical flare lagging the X-ray maximum by ~ 10 days. 
We also notice that the fractional variability of the source is com- 
parable at soft (22.5 %) and hard (25 %) X-rays. Interestingly, 
we do not find any significant correlation among the X-ray spec- 
tral index and flux. This may be due to the poor data sampling 
or may be intrinsic to the source. The concave shape of the X- 
ray spectrum (see Section 13.5) , suggests that the X-ray emis- 
sion shows a combination of synchrotron and inverse-Compton 
mechanisms which could prevent the source from exhibiting any 
steepening or hardening trend during the flare. 

A simil ar orphan X - ray flar e was also observed in the blazar 
3C 279 by lAbdo et al.l (l2010bl) with X-ray flaring activity lag- 
ging optical - GeV flares by 60 days. The authors argued that 
X-ray photons are produ ced further down to the jet compared 
to optical - GeV photons. iHavashida et al.l (1201 2h argued in the 
context of a two component model; the X-ray flare is produced 
by the low-frequency component which is less variable com- 
pared to the high-frequency component. Although we do not 
completely understand the origin of the orphan X-ray flare in 
S5 0716+714, it appears possible that the X-ray emission is not 
co-spatial with the optical/y-ray emission in this event. We no- 
tice some low level flux activity (mini flare, say R7) in between 
the two major radio flares ("R6" and "R8", see Fig. |2]i. While 
modeling the radio spectra of the source, we also noticed that a 
two-component model better describes the synchrotron spectra 
of the source over this period. This indicates that either mul- 
tiple shocks are hitting the emission region which at first pro- 
duces the major flare "05/G5-R6", then "X6/06-R7" and later 
"07/G7-R8", or the radiation is contributed by two synchrotron 
components with the low-frequency component producing the 
X-ray flare. 



5. Conclusions 

In this paper we presented the results of the radio to y-ray mon- 
itoring of S5 0716+714 from April 2007 to January 2011. The 
source was very active at optical and higher frequencies. Two 
major radio outbursts were observed during this high activity pe- 
riod. From the rapid rise and decay, we derive variability bright- 
ness temperatures exceeding the IC limit, which at least for mm 
flares is a very unique behavior 

A long-term variability trend (~ 350 days timescale) is visi- 
ble in the optical light curves which is superimposed with repet- 
itive variations on shorter time scales (~ 60 day). A comparison 
of the various flaring episodes of S5 0716+714 strongly indi- 
cates a one-to-one correlation between the strength of the y-ray 
emission and the strength of the optical emission. A quadratic 
dependence of the amplitude of the y-ray variabihty with respect 
to that of the optical favors an SSC explanation. 

The high-energy (optical - GeV) flares propagate down to 
radio frequencies with a time delay of ~ 65 days following a 
power-law dependence on frequency with a slope ~ 0.3. This 
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indicates that opacity plays a key role in producing time delays 
among light curves at optically thin and thick wavelengths. Since 
the radio outbursts are low-peaking flares, such a long time lag is 
only possible in the case of optical flares being the precursors of 
radio ones. The evolution of the radio flares are in agreement 
with t he generalized shock model proposed by IValtaoia et al.l 
(Il992h . The evolution of the flare in the turnover frequency - 
turnover flux density ( v,„ - S ,„) plane shows a very steep rise and 
decay over the Compton and adiabatic stages with a slope too 
high to be expected from intrinsic variations, requiring an addi- 
tional Doppler factor variation along the jet. For the two flares, 
we notice that 6 changes as R^-^ during the rise and as ^ during 
the decay of R6 flare. The evolution of the R8 flare is governed 
by 5 oc "* during the rising phase and 6 oc R^^ ^ during the 
decay phase of the flare. 

An orphan X-ray flare is observed ~ 50 days after the ma- 
jor optical - GeV flares. The detection of an isolated X-ray 
flare challenges the simple one-zone emission models, render- 
ing them too simple. The lack of substantial observations over 
the flaring epoch makes it even more complicated to understand. 
We found that this flare has equal contributions from both the 
synchrotron and the high-energy (inverse-Compton in a leptonic 
model interpretation) emission mechanisms. 

We model the broadband SEDs of the source using two dif- 
ferent versions of leptonic models: a pure SSC and SSC+EC. We 
found that the low activity states of the source are well described 
by a pure SSC model while an EC contribution is required to re- 
produce the SEDs for high states. The SSC+EC model returns 
magnetic field parameter value closer to equipartition, providing 
a satisfactory description of the broadband SEDs. We found that 
satisfactory model fits can be achieved if the external radiation 
field is dominated by Ly-a emission from the broad-line region. 
This model nicely describes the broadband SEDs of the source 
at optical and higher frequencies, but under-predicts the cm - 
mm spectra at least for few time periods. A separate synchrotron 
component seems required to fit the cm - mm radio fluxes. This 
may also provide a hint towards the origin of the orphan X-ray 
flare. 
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